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Chapter  1 

INTRODUCTION 


1.1  Overview  of  the  simulation  code  HPHall-2 

HPHall-2[3]  is  a  two-dimensional(2D),  non-stationary,  hybrid  (PIC-fluid)  and  quasineutral  code 
for  Hall  Effect  Thruster  (HET)  simulation.  Heavy  species(ions  and  neutrals)  are  modelled  with 
a  particle-in-cell  plus  Monte-Carlo  Collisions(PIC-MCC)  method  while  electrons  are  modelled 
macroscopically.  This  procedure  is  a  good  trade-off  between  fully  fluid  models,  which  are  not 
well  suited  for  Hall  thrusters  where  particle  densities  are  too  low,  and  fully  kinetic  models  where 
computational  time  is  enormous.  In  addition,  HPHall-2  assumes  cylindrical  symmetry  and  lat¬ 
eral  ceramic  walls  and  uses  an  externally-imposed  magnetic  held,  B.  Plasma  quasi- neutrality  is 
based  on  the  Debye  length  being  much  smaller  than  typical  lengths  of  Hall  thrusters.  Therefore, 
electron  density  and  singly  and  doubly  charged  ion  densities  are  directly  related.  However,  due 
to  boundary  conditions  the  quasi-neutrality  assumption  is  not  valid  in  very  thin  layers,  called 
sheaths,  tied  to  the  walls.  As  a  consequence  the  boundaries  of  the  simulation  domain  are  not 
the  walls  but  the  transition  surfaces  to  the  space-charge  sheaths  tied  to  the  chamber  walls. 
This  requires  to  solve  independently  in  a  sheath  subcode  the  different  sheath  problems  in  order 
to  determine  the  correct  conditions  at  the  boundaries  of  the  quasineutral  domain.  Figure  1.1 
shows  the  different  elements  of  the  simulated  thruster,  the  mesh  used  by  the  PIC  part  of  the 
code,  the  magnetic  held  profile  at  the  channel  median,  and  the  magnetic  streamsurfaces  used 
by  the  electron  subcode. 

PIC  methods  simulate  species  as  a  group  of  macroparticles,  each  one  of  them  representing 
a  large  number  of  particles.  Macro-particle  motion  is  computed  by  integrating  Newton's  law 
while  collision  processes,  such  as  ionization  or  charge  exchange  collisions,  are  modelled  using 
a  Monte-Carlo  Collision  algorithm.  Forces  acting  on  particles  are  due  to  the  electric  field,  the 
externally  applied  magnetic  held  and  ion-neutral  collisions.  Macroscopic  magnitudes  for  ions 
and  neutrals  (density,  flux,  temperature,...)  are  computed  by  weighting  particle  magnitudes 
at  every  node  of  the  mesh.  Typical  numbers  are  macro-particles  30.000  per  species  (neutrals, 
singly  charged  ions  and  doubly  charged  ions)  and  about  1000  cells  in  the  mesh. 

The  strong  magnetic  held  causes  a  high  anisotropy  in  the  electron  macroscopic  dynamics, 
which  allows  to  treat  the  directions  parallel  and  perpendicular  to  B  in  different  manner.  In 
particular,  Boltzmann  equilibrium  is  assumed  along  magnetic  lines  while  an  Ohm’s  law  (with 
pressure  effects)  is  used  in  the  perpendicular  direction.  The  2D  energy  equation  is  transformed 
into  a  ID  problem  by  first  integrating  electron  equations  along  the  magnetic  held  lines. 
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Figure  1.1;  Elements  of  the  simulated  SPT-100  type  simulated  thruster  and  PlC-mesh(a),  reference 
magnetic  field  at  the  chamber  median,  and  are  the  radial  and  axial  components  respectively(b) 
and  magnetic  streamsurfaces  used  in  the  electron-mesh(c). 


The  main  outputs  of  the  PIC  sub-code  are  densities  and  particle  fluxes  computed  by  weight¬ 
ing  algorithms.  These  values  are  needed  by  the  electron  sub-code  that,  in  turn,  computes  the 
electron  temperature,  necessary  to  evaluate  collisional  processes,  and  the  electric  potential, 
which  dehnes  the  electric  held  to  be  applied  on  ions. 

1.2  Objectives  of  the  project 

This  project  is  aimed  to  improve  HPHall-2  (as  delivered  in  grant  FA8655-04-1-3003)  in  two 
categories  of  objectives: 

1.  A  better  modelization  and  more  accurate  algorithms  of  the  main  phenomena  of  the  dis¬ 
charge.  This  includes: 

(a)  Detailed  treatment  of  the  near  anode  region  for  oblique  magnetic  helds.  The  diffi¬ 
culty  lies  in  matching  the  ID  electron  equations  at  each  magnetic  line  to  the  anode 
sheath  boundary,  which  is  not  a  magnetic  streamline.  The  handling  of  boundary 
conditions  for  the  case  of  an  anode  sheath,  that  can  be  electron-repelling  in  certain 
zones  and  ion-repelling  in  others,  is  an  important  issue,  partially  unsolved  yet. 

(b)  Implementation  of  more  accurate  sheath  models. 

(c)  Improvement  of  algorithms  for  moving  and  injecting  particles  in  the  PIC  subcode. 
The  objective  is  a  reduction  of  temporal  errors  by  one-order  of  magnitude.  Apart 
from  obtaining  a  higher  accuracy,  this  will  allow  to  carry  out  long  runs  and  numerically- 
reliable  parametric  studies. 

(d)  Improvement  of  magnitude  weighting  at  the  domain  boundaries. 

(e)  The  development  of  a  more  accurate  electron  subcode.  A  new  differential  formulation 
of  the  electron  equations,  totally  based  on  curvilinear  magnetic  coordinates  is  being 
derived.  Modihed  algorithms  try  to  increase  the  spatial  and  temporal  accuracy  of 
the  solution.  Effects  previously  neglected  like  electric  work  parallel  to  B-lines  are 
being  included.  A  modihcation  of  the  constant-temperature,  Boltzmann  equilibrium 
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model  in  the  B-parallel  direction  and  a  more  detailed  model  of  anomalous  diffusion 
are  pending  on  advances  on  parallel  researches. 

2.  To  increase  the  capabilities  of  the  code. 

(a)  High  specihc  impulse  performances.  Double  ions  become  a  signihcant  population  at 
high  discharge  voltages.  The  new  code  includes  double  ions  at  a  separate  species. 
Effects  of  double  ions  on  the  sheath  conditions  are  included  too.  A  higher  flexibility 
on  the  magnetic  meshes  to  be  managed  by  the  electron  subcode  is  necessary. 

(b)  Two-stage  discharges.  This  consists  mainly  on  considering  the  presence  of  internal 
electrodes  in  certain  regions  of  the  inner  or  outer  walls.  The  implications  for  the 
code  are:  (1)  another  major  modihcation  on  the  electron  subcode,  since  the  dis¬ 
charge  current  cannot  be  considered  constant  along  all  magnetic  streamlines;  (2) 
the  formulation  of  a  sheath  model  for  different  electrode  types,  in  particular  for 
highly-emissive  ones. 


1.3  Organization  of  the  report 

The  changes  implemented  in  the  code  are  detailed  in  chapters  2  to  6.  Chapter  2  is  devoted 
to  changes  in  the  general  structure  of  the  code  and  in  the  pre-process  subroutines  (mesh  and 
wall  dehnitions,  etcetera).  Chapters  3  to  5  are  devoted  to  explain  the  theoretical  advances 
that  have  been  implemented  in  the  sheath,  PIC,  and  electron  subcodes,  respectively.  Chapter 
6  illustrates  some  of  the  capabilities  of  the  new  version.  Finally  Annexes  A  and  B  present 
on-going  theoretical  advances  that  have  not  been  implemented  in  the  code  yet. 
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Chapter  2 

GENERAL  FEATURES  OF  THE  CODE 


2.1  Improved  code  structure 

HPHall-2  has  been  modified  in  order  to  make  it  easy-to-use  by  people  not  familiar  with  the 
computational  code  itself.  It  has  been  divided  into  several  applications  each  one  of  them 
designed  to  perform  a  different  task.  In  this  manner,  there  are  separate  programs  to: 


-  load  a  spatial  grid  computed  with  TecPlot  and  generate  a  file  in  the  format  required  by 
HPHall. 

-  load  a  magnetic  field  computed  with  Maxwell  and  generate  a  file  in  the  format  required 
by  HPHall. 

-  generate  a  magnetic  field  file  in  the  format  required  by  HPHall  based  on  a  theoretical  law 
that  can  be  defined  by  the  user. 

-  generate  files  containing  tables  that  represent  the  wall  interaction  models. 

-  execute  HPHall  using  the  different  files  generated  in  the  previous  steps  as  well  as  other 
physical  and  numerical  parameters  contained  in  a  unique  file. 


Apart  from  the  files  generated  by  the  previous  applications,  HPHall  needs  one  input  file 
where  both  physical  constants  (e.g.  Xenon  properties)  and  numerical  parameters  (e.g.  ion  and 
electron  time  steps)  are  defined.  In  this  way,  the  whole  execution  can  be  set  up  from  one  single 
file.  The  main  variables  that  must  be  defined  in  this  file  are  shown  in  several  tables  at  the  end 
of  this  chapter. 

HPHall  generates  as  output  self-explanatory  data  files  that  can  be  analyzed  either  using 
TecPlot  or  a  simple  text  editor.  For  instance,  the  program  generates  an  output  file  containing 
the  main  parameters  regarding  thruster  operation;  efficiency,  thrust,  power,  discharge  voltage, 
etc.  Besides,  several  files  are  generated  containing  different  spatial  and  temporal  profiles  for 
both  wall  properties  and  plasma  magnitudes  inside  the  thruster.  All  this  files  are  generated  in 
an  independent  directory  so  that  results  from  different  executions  can  be  compared  easily. 
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2.2  New  algorithm  to  compute  the  magnetic  grid 

HPHall-2  solves  electron  equations  using  a  ’magnetic  grid’  derived  from  the  magnetic  field 
topology  and  the  spatial  grid.  However,  it  fails  to  compute  properly  this  magnetic  grid  and 
consequently,  errors  are  propagated  to  the  electron  equations. 

In  order  to  solve  that  problem  a  new  algorithm  has  been  implemented  that  ensures  the 
magnetic  grid  is  computed  consistently.  In  fact,  this  new  algorithm  is  more  consistent  with  the 
geometrical  grid  because  it  reproduces  any  grid  rehnement  that  may  appear  near  the  lateral 
walls. 

Figure(2.1)  depicts  a  comparison  between  the  magnetic  grids  obtained  using  the  algorithm 
in  HPHall-2  and  the  new  one.  HPHall-2  fails  to  reproduce  properly  the  lateral  walls  as  well 
as  the  near  plume  region.  This  fact  causes  wall  properties  are  not  evaluated  in  the  wall  but 
inside  the  computational  domain.  Thus  plasma-wall  interaction  properties  are  not  computed 
accurately. 


Figure  2.1:  Spatial(left)  and  magnetic(right)  grids  computed  by  HPHall-2  (upper  figures)  for  the 
SPT-70  thruster  and  computed  by  the  new  version  with  the  new  algorithm  (lower  figures)  for  the 
SPT-100  thruster. 

In  the  upper  hgures  (i.e.,  SPT-70  hgures)  it  can  be  observed  that  the  algorithm  implemented 
in  HPHall-2  does  not  compute  accurately  the  simulation  domain  and  the  near  plume  region 
upper  limits.  In  fact,  the  upper  lateral  wall  should  be  at  0.035  m  from  the  centerline  (see 
geometric  grid)  while  the  magnetic  grid  does  not  reach  that  value  at  any  point.  On  the  other 
hand,  the  new  algorithm  (i.e.,  SPT-100  hgures)  generates  a  magnetic  grid  completely  equivalent 
to  the  geometric  grid. 
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2.3  Channel  geometry  with  curved  walls 

Another  interesting  improvement  is  the  added  capability  of  the  new  code  to  deal  with  curved 
walls.  This  is  interesting  in  case  the  influence  of  wall  erosion  is  to  be  analyzed  on  thruster 
performance. 

Figure  2.2  shows  the  results  obtained  with  the  new  code  for  both  a  curved  and  non-curved 
wall.  Obviously,  to  make  full  profit  of  this  improvement  a  sputtering  algorithm  must  be  im¬ 
plemented  to  compute  the  wall  erosion  rate.  However,  this  is  out  of  the  scope  of  the  current 
project. 

As  expected,  as  the  cross  section  area  increases  near  the  exit,  the  plasma  density  decreases, 
the  ion  flux  is  spread  and  the  temperature  decreases. 


Figure  2.2:  2D  maps  comparing  results  for  curved  and  non-curved  walls. 


2.4  Local  definition  of  wall  properties 

In  order  to  increase  the  capabilities  of  the  code  an  additional  feature  has  been  added.  Now  it  is 
possible  to  dehne  as  many  different  materials  (or  even  different  models  for  the  same  material) 
as  desired.  This  is  controlled  in  HPHall  main  input  hie.  In  this  hie  there  is  a  group  of  variables 
that  can  be  used  to  dehne  the  diherent  material  geometries,  properties,  models,  etc.  At  the  end 
of  the  chapter  the  main  values  to  be  dehned  are  presented  in  a  table.  The  diherent  wall  segments 
made  of  a  given  material  must  be  dehned  with  r  and  z  coordinates.  Besides,  each  material  must 
have  specihed  the  material  name,  type  and  hie  containing  the  plasma-wall  interaction  table. 

In  the  last  chapter  results  are  presented  for  a  two-stage  thruster  where  an  intermediate 
electrode  (modeled  as  an  ideal  electron  emitter)  is  presented.  This  is  a  good  example  of  the 
hexibility  allowed  with  this  new  feature. 
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Parameter  name 

Value 

Description 

B_FILE 

mag_SPT100.dat 

magnetic  field  input  file 

GRID_TECPLOT_FILE 

grid_SPT100.dat 

grid  input  file 

Table  2.1;  HPHall  input  parameters  (INPUT  FILES). 


Parameter  name 

Value 

Description 

RUN_MODE 

0 

0  =  normal,  1  =  neut.,  2  =  ions,  3  =  particles 

N_ITS 

10000 

number  of  iterations 

TRANSIENT_ITS 

0 

number  of  iterations  considered  as  transient 

SAVE_ITS 

2000 

save  data  every  x  iterations 

CONCATENATE 

0 

concatenate  averaged  values 

N_CONCAT_ITS 

0 

number  of  previous  iterations  to  concatenate 

Table  2.2;  HPHall  input  parameters  (PROGRAM  FLOW). 


Parameter  name 

Value 

Description 

MOVIES 

0 

save  movies 

MOVIE_ITS 

50 

save  movie  frame  every  x  iterations 

MOVIE_FLUID 

1 

save  fluid  movie 

MOVIE_PARTICLE 

0 

save  particle  movie 

N_IONS_MOVIE 

1000 

number  of  ions  for  the  particle  movie 

N_NEUTRALS_MOVIE 

1000 

number  of  neutrals  for  the  particle  movie 

Table  2.3;  HPHall  input  parameters  (MOVIES). 


Parameter  name 

Value 

Description 

SAVE_ELECTRON_INFO 

0 

save  electron  information  (y/n) 

SAVE_ELECTRON_ITS 

1 

save  electron  information  every  x  electron  iterations 

Table  2.4;  HPHall  input  parameters  (ELECTRON  MOVIE  AND  MAP). 
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Parameter  name 

Typical  value 

Description 

PROPELLANT 

2 

ARGON;  1  XENON;  2 

MDOT 

4.8e-6 

mass  flow  rate  in  kg/s 

PPU_TYPE 

1 

GONSTANT_VOLTAGE=l 
GONSTANT  GURRENT=2 
GONSTANT_POWER=3 

V_DISCHARGE 

300.0 

discharge  voltage  (for  constant  voltage) 

I_DISGHARGE 

6.5 

total  current  at  cathode  (for  constant  current) 

P_DISGHARGE 

2000 

total  power  (for  constant  power) 

Table  2.5;  HPHall  input  parameters  (OPERATION  PARAMETERS). 


Parameter  name 

Value 

Description 

WALL_FLUX_V_BBG 

1 

wall  normal  velocity  is  dehned  using  the  wall  flux 

GORREGTED_WEIGHTING 

0 

new  algorithms  to  weight  at  the  boundaries 

SURFAGE_WEIGHTING 

1 

new  algorithms  to  weight  at  boundaries 

BOHMFORGING 

0 

apply  Bohm  forcing  algorithm 

ANODE_BOHM_FORGING 

0 

anode  is  included  in  bohm  forcing 

NEW_IONIZATION_LAW 

1 

new  ionization  law  for  particle  number  control 

NEW_GATHODE 

0 

new  cathode  model 

Table  2.6:  HPHall  input  parameters  (MODELLING  OPTIONS). 
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Parameter  name 

Value 

Description 

N_MATS 

3 

number  of  different  materials 

WALL_TYPE_1 

1 

anode(l),  ceramic(2),  electrode(3),  cathode(4) 

MATERIAL_NAME_1 

Anode  ERF  200eV 

material-model  name 

SEE_C0EFF_P_1 

0.000 

SEE  law  coefficient  p  (meaningless  for  anode) 

SEE_E1_1 

0.000 

SEE  law  energy  (J)  (meaningless  for  anode) 

N_SEGMENTS_1 

1 

segment  number  for  material  1 

Z_BEGIN_1_1 

0.0000 

z-coordinate  of  segment  1  starting  point 

R_BEGIN_1_1 

0.0350 

r-coordinate  of  segment  1  starting  point 

Z_END_1_1 

0.0000 

z-coordinate  of  segment  1  ending  point 

R_END_1_1 

0.0500 

r-coordinate  of  segment  1  ending  point 

P0TENTIAL_1_1 

0.0 

potential(V)  (meaningless  for  ceramic  material) 

WALL_TYPE_2 

2 

MATERIAL_NAME_2 

BN_ERF_200eV 

SEE_G0EFF_P_2 

0.576 

SEE_E1_2 

7,63E-15 

N_SEGMENTS_2 

4 

Z_BEGIN_1_2 

0.0000 

R_BEGIN_1_2 

0.035 

Z_END_1_2 

0.0250 

R_END_1_2 

0.035 

P0TENTIAL_1_2 

0.0 

Z_BEGIN_2_2 

0.0000 

R_BEGIN_2_2 

0.050 

Z_END_2_2 

0.025 

R_END_2_2 

0.050 

POTENTIAL_2_2 

0.0 

Z_BEGIN_3_2 

0.025 

R_BEGIN_3_2 

0.035 

Z_END_3_2 

0.025 

R_END_3_2 

0.00977 

POTENTIAL_3_2 

0.0 

Z_BEGIN_4_2 

0.025 

R_BEGIN_4_2 

0.050 

Z_END_4_2 

0.025 

R_END_4_2 

0.0693 

POTENTIAL_4_2 

0.0 

WALL_TYPE_3 

3 

MATERIAL_NAME_3 

electrode 

SEE_G0EFF_P_3 

0.576 

SEE_E1_3 

7,63E-15 

N_SEGMENTS_3 

0 

Z_BEGIN_1_3 

0.0250 

R_BEGIN_1_3 

0.0580 

Z_END_1_3 

0.0250 

R_END_1_3 

0.0600 

P0TENTIAL_1_3 

0.0 

Table  2.7;  HPHall  material  definition 


Chapter  3 

IMPROVEMENTS  IN  SHEATH 
SUBCODE 


HPHall  is  a  quasineutral  code,  the  simulation  boundaries  are  not  the  walls  but  the  transition 
surfaces  to  the  space-charge  sheaths  tied  to  the  chamber  walls.  This  requires  to  solve  indepen¬ 
dently  the  different  sheath  problems  in  order  to  dehne  the  correct  conditions  at  the  boundaries 
of  the  quasineutral  domain.  This  chapter  explains  the  different  models  used  for  the  differ¬ 
ent  kinds  of  sheaths  that  can  appear  inside  a  Hall  thruster;  anode,  ceramic-wall  sheaths  and 
intermediate  electrodes. 


3.1  Anode  sheath 

A  simple  sketch  of  the  anode  is  represented  in  hgure  3.1(a).  Because  of  its  thinness,  of  the 
order  of  the  Debye  length,  the  sheath  is  considered  collisionless,  quasisteady,  and,  in  principle, 
negative  (meaning  electron- repelling).  Since  the  wall  is  considered  metallic,  the  role  of  this 
sheath  is  to  adjust  the  small,  diffusive  electron  flux  coming  from  the  quasineutral  plasma  to 
the  electron  flux  collected  by  the  (perfectly-absorbing)  anode.  Hence,  the  main  sheath  charac¬ 
teristics  depend  only  on  electrons. 

The  anode  sheath  model  used  by  Fife,  Ahedo  et  ah  [4,  5]  and  other  researchers  assumes  a 
full  Maxwellian  distribution  of  electrons  at  the  sheath  entrance  (point  B).  This  model  is  inac¬ 
curate  since  the  anode  does  not  return  back  the  electrons  reaching  it.  A  drifted  Maxwellian,  as 
assumed  by  Ahedo [6],  is  incorrect  too,  since  the  distribution  must  be  symmetric  on  the  normal 
velocity  for  electrons  that  are  reflected  back  within  the  sheath. 


The  correct  choice,  within  the  ’family’  of  Maxwellian  distributions,  is  a  truncated,  zero- 
drift  distribution.  Thus,  taking  into  account  the  energy  conservation  of  the  electrons  within 
the  sheath,  one  has 


/e(u±,U||,0)  =  n* 
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27rT* 


(3/2) 

exp 


exp 


me{v\  +  vh 


2U 


H[v^  + 


'2e(0  -  (Pa) 


TTle 


(3,1) 


23 


24 


CHAPTER  3.  IMPROVEMENTS  IN  SHEATH  SUBCODE 


Plasma- sheath 

Anode  transition  x 

A  B 


.  gx 


Non 

quasineutral  Quasineutral 
sheath  plasma 


Figure  3.1;  (a)  Sketch  of  the  anode  sheath,  (b)  The  truncated  Maxwellian  distribution  at  the  sheath 
side  of  B,  and  the  ’equivalent’  drifted  Maxwellian  distribution  corresponding  to  macroscopic  magnitudes 
at  the  quasineutral  side  of  B.  In  these  figures  the  x  axis  represents  the  direction  perpendicular  to  the 
wall 


Here,  fe  is  the  electron  distribution  function  for  the  velocity  component  perpendicular  to  the 
anode  surface,  v±,  the  parallel  velocity  component,  ny,  and  the  electric  potential,  0,  which  plays 
the  role  of  the  spatial  variable;  H{v)  is  the  Heaviside  step  function  of  n;  and  T*  and  n*  are 
constants  characterizing  the  distribution.  These  constants  are  the  temperature  and  density  at 
B  only  in  the  limiting  case  of  a  non-truncated  distribution. 


Let  us  consider  a  one-dimensional  problem  where  plasma  magnitudes  are  constant  along 
the  whole  surface  transition  B.  If  (pAB  is  the  sheath  potential  fall  and  taking  the  appropriate 
moments  of  fe,  the  density,  particle  flux,  temperature,  and  total  energy  flux  to  the  sheath  at 
point  B  are 


l  +  eiii^yecpAB/T^) 

rieB  =  - X - 


(3.2) 


g±eB  — 


e(j)AB\ 


(3.3) 


I -  -  f  _  e4>AB\ 

TeB  .  ,  _ V  T,  ) 

r.  ST.  T.  V  T.  A+erf(yeWr.)’ 

Q±eB  —  gXeBif^Ti,  -f-  C^ab),  (3.5) 

respectively,  with  u±e  =  Qxe/ne  the  macroscopic  electron  velocity  and  erf(a;)  the  error  function. 


These  magnitudes  must  match  with  those  known  from  the  quasineutral  macroscopic  model 
at  the  sheath  transition.  Therefore,  these  four  equations  determine  the  potential  fall  (pAB,  the 
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Figure  3.2:  Comparison  of  anode  models  in  terms  of  the  non-dimensional  electron  density  (a),  non- 
dimensional  electron  temperature  (b),  electron  Mach  number  (c),  non-dimensional  electron  flux  to 
anode  (d),  non-dimensional  total  energy  flux  to  the  anode  (e)  and  non-dimensional  heat  flux  to  the 
anode  (f). 


constants  of  the  distribution  function,  n*,  T*,  and  the  heat  conduction  flux  to  the  sheath; 

(l±eB  =  ~  ■^9±eBTeB 

Figure  3.2  compares  the  sheath  parameters  for  the  new  and  old  (i.e.  full  Maxwellian) 
sheath  models.  As  expected,  both  models  coincide  exactly  for  an  inhnite  sheath  potential  fall 
and  approximately  for  e(j)AB/TeB  >  3.  The  old  model [7,  8]  placed  the  transition  to  a  no-sheath 
regime  (i.e.  4>ab=^)  at  M^b  =  (27r)“^/^  ~  0.4,  with  Me^  =  \u±eB\/ fhe  electron 
Mach  number.  The  present  model  places  it  at  M^b  ~  0.9,  much  closer  to  M^b  =  1  where 
the  quasineutral,  macroscopic  model  of  Ahedo  [6]  places  the  transition  from  a  negative-  to  a 
positive-  anode  sheath  (for  near-cold  ions).  [A  perfect  parametric  matching  of  the  quasineutral 
and  kinetic  models  is  unlikely  since  they  have  been  postulated  independently.]  Notice  that  the 
increment  of  S'j.e/'n-eB  and  Q^Ib/'^cB  as  the  sheath  vanishes  is  due  to  the  increase  of  u±eB  and 
decrease  of  n^B)  and  not  to  changes  of  g±e-  Finally,  another  merit  of  the  new  model  is  that  heat 
conduction  to  the  anode  is  positive  always;  the  old  model  predicted  a  negative  heat  conduction 
to  the  anode  in  the  no-sheath  limit,  which  did  not  look  satisfactory. 

In  this  analysis  ions  have  played  no  role.  The  Bohm  condition  states  that  the  ion  flux  must 
enter  sonic-supersonically  into  a  negative  sheath.  The  fulhllment  of  this  condition  pertains  to 
the  PIC  subcode  and,  particularly,  to  the  topic  of  plasma  weighting  at  the  boundaries. 

A  recent  improvement  of  this  anode-sheath  model,  which  places  the  vanishing  of  the  negative 
sheath  exactly  at  M^b  =  1  is  presented  in  Annex  A,  but  has  not  been  implemented  in  the  code 
yet. 
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3.2  Ceramic  wall  sheath 

The  sheath  around  ceramic  walls  is  more  complex,  since  (i)  both  the  electron  and  ion  dynamics 
must  be  taken  into  account,  (ii)  there  is  secondary  electron  emission  (SEE)  by  electron  impact 
in  the  walls,  and  (hi)  space-charge  saturation  is  possible [9]. 

Improvements  on  the  sheath  model  are  being  implemented  on  three  points.  First,  the  use  of 
truncated  Maxwellian  distribution  for  electrons,  instead  of  the  full  Maxwellian  distribution  used 
originally  by  Hobbs  and  Wesson[9]  and  adopted  by  Ahedo[10].  Second,  the  inclusion  of  doubly- 
charged  ions  as  an  additional  species,  since  they  constitute  a  non-marginal  population  for  high 
specihc  impulse  operation  [11,  1],  And  third,  the  consideration  of  a  supersonic  ion  flux  as  a 
likely  sheath  entrance  condition,  because  of  the  discharge  fluctuations.  With  respect  to  this  last 
point,  we  know  that  the  Bohm  condition  for  the  sheath  states  that  ions  must  enter  sonically  or 
supersonically.  At  the  same  time,  the  stationary  (and  macroscopic)  solution  for  the  quasineutral 
plasma  states  that  ions  must  enter  into  the  sheath  subsonically  or  sonically.  However,  supersonic 
ion  fluxes  are  possible  in  a  time-fluctuating  quasineutral  solution,  as  Parra  et  al.[12]  showed. 
Finally  observe  that  the  sheath  remains  quasi-steady  for  a  fluctuating  discharge. 
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quasineutral 

sheath 


Quasineutral 

plasma 


W  Wall 


Sri+  gri++  Srs  Plasma-sheath 
^  transition 

§ri+  S  ri-i-4-  g  re  1 
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Figure  3.3:  Sketch  of  ceramic  wall  fluxes  and  nomenclature. 

The  fundamentals  of  the  improved  model  are  those  set  by  Ahedo[10]  for  total  thermalization 
except,  of  course,  for  the  aspects  we  want  to  improve.  The  issue  of  partial  thermalization  of 
the  electron  distribution  is  not  treated  here,  thus,  total  thermalization  is  assumed.  We  omit 
here  the  derivation  of  the  model  and  just  include  some  notes  on  it.  A  sketch  of  the  ceramic  wall 
sheath  is  shown  in  hgure  3.3.  The  plasma  species  within  the  sheath  are:  primary  electrons(p) 
coming  from  the  quasineutral  plasma,  secondary  electrons(s)  emitted  by  the  wall,  singly-charged 
ions(f+)  and  doubly-charged  ions(f++).  Primary  and  secondary  electrons  match  into  a  single 
electron  population  (e)  outside  the  sheath.  Quasi  cold  beams  are  assumed  for  the  non-conhned 
populations  {i+,  i++,  and  s).  The  zero  current  condition  leads  to 

9ri+W  T  ‘^9ri-\ — \-W  (f  ^w)9rpWi  (^’S) 

where  Qr  mean  particle  flux  perpendicular  to  the  (radial)  wall,  and  5^  is  the  effective  SEE  yield. 

Plasma  quasineutrality  at  the  sheath  entrance  yields 

^sQ  T  ^pQ  —  ^eQ  ^i+Q  T  ^^i++Q  (3-7) 

and  the  average  ion  charge  number  for  ions  at  point  Q  is 

^iQ  -  -  -  - ^ - • 

^iQ  '^i+Q  T  '^i++Q 


(3.8) 
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Figure  3.4:  Ceramic  material  sheath  model  non-dimensional  results  and  influence  of  the  supersonic 
ion  flow  on  the  sheath  properties. 


When  there  is  only  one  ion  beam,  the  Bohm  condition  of  the  sheath  is  explicit  for  the  ion 
macroscopic  velocity,  Eq.(18)  of  Ref.  [10].  In  the  presence  of  singly-charged  and  doubly-charged 
ion  beams  the  Bohm  condition  yields  only  one  relation  between  the  two  beam  velocities.  A 
correct  closure  of  the  problem  would  require  to  involve  the  quasineutral  solution,  and  this  would 
be  very  costly.  Two  simple,  alternative  procedures  are  being  studied.  First,  to  take 

^ri++Q  (2’^) 

as  if  the  two  ion  beams  were  accelerated  freely  in  the  presheath,  with  no  ion  production.  Second, 
to  consider  a  single  ion  beam  of  charge  ^iQ-  Results  here  are  presented  for  the  hrst  alternative. 
The  sheath  potential  fall,  0vvq,  turns  out  to  verify 


^4>wq 

T* 


In 


-F  ln(l  —  -\-  In 


In  (2^2  -  1  -  2(^2  -  -  In  M^q 


In 


1  -h  erf  ( ^Je(t)wQ/T^ 


(3.10) 


where  the  second  line  groups  the  effects  newly  introduced:  doubly-charge  ions,  ZiQ  >  1,  su¬ 
personic  ion  flux  MiQ  =  u j_|_Q /  a/ TeQ / m j ,  and  truncated  Maxwellian  distribution  (based  on 
’temperature’  T*)  for  primary  electrons  (yielding  the  term  with  the  error  function). 

This  new  model  depends  on  three  dimensionless  parameters  coming  from  the  quasineutral 
solution  (i.e.  from  the  simulation  of  the  2D  discharge)  and  the  wall  material:  Teq/Ti  (with 
Ti  the  temperature  for  a  100%  SEE  yield[10])  Mjq,  and  Ziq.  Figure  (3.4)  shows  that  the  ion 
Mach  number  is  found  to  affect  signihcantly  the  potential  fall  and,  thus  the  energies  of  ions  and 
electrons  impacting  the  wall.  The  influence  of  the  ion  charge  number  can  be  interpreted  as  an 
increase  of  ion  Mach  number,  as  Eq.  (3.10)  shows.  Since,  according  to  Ref.[l],  Ziq  is  between 
1  and  1.1  its  effect  is  small. 

In  the  present  model  the  charge  saturation  limit (CSL)  is  not  universal,  but  depends  on  two 
parameters,  Miq  and  Z^q.  This  is  illustrated  in  hgures  (3.5)  and  (3.6).  Both  new  effects  tend 
to  reduce  the  value  of  the  SEE  yield  at  the  CSL,  thus  reducing  the  electron  energy  losses  to 
the  walls. 

Finally,  the  use  of  a  truncated  distribution  function  for  primary  electrons  at  Q  introduces 
modihcations  only  at  the  CSL  or  close  to  it.  For  the  basic  case,  M^q  =  1  and  Z^q  =  1, 
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Figure  3.5;  Influence  of  the  ion  Mach  number  on  the  properties  at  the  saturation  charge  limit.  SEE 
yield  at  CSL(a),  non-dimensional  temperature(b),  non-dimensional  sheath  potential  fall(c)  and  non- 
dimensional  primary-electron  density  (d). 


Figure  3.6;  Influence  of  the  equivalent  charge  number  on  the  charge  saturation  limit. 


Hobbs  and  Wesson  had  =0.983,  whereas  now  we  hud  =0.9865.  One  can  be  tempted 
to  neglect  this  0.3%  difference  on  unless  one  realizes  that  the  important  magnitude  is 

1/(1  —  5^^^)  which  measures  particle  and  energy  fluxes  to  the  walls.  Then,  one  hnds  that  the 
differences  in  energy  deposition  with  Hobbs- Wesson  amount  to  a  25%  (larger  in  the  present 
case) . 

Figure  3.7  displays  a  comparison  of  the  results  obtained  with  the  simulation  code  of  Ref.[l] 
using  the  old  and  new  models  for  the  lateral  sheaths. 
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Figure  3.7;  Plasma  characteristics  at  the  outer  wall  with  the  new  (solid  lines)  and  old  (dashed  lines) 
models  for  the  lateral  sheaths.  The  simulation  data  is  the  same  used  in  Ref.[l]. 


3.3.  EFFECTS  OF  ELECTRON  PARTIAL  THERMALIZATION 


29 


3.3  Effects  of  electron  partial  thermalization 

The  sheath  model  at  the  lateral  walls  of  the  previous  section  is  based  on  the  assumption  that 
the  thermalization  mean  free-path  for  electrons  is  smaller  than  the  channel  width.  In  this  case, 
the  emitted  beams  of  secondary  electrons  thermalize  with  the  bulk  electron  population  and  the 
tail  of  electrons  collected  at  one  wall  is  quickly  replenished.  As  a  consequence,  the  velocity 
distribution  function  of  ’primary’  electrons  flowing  into  the  sheath  edge  can  be  considered 
semi-Maxwellian. 

In  practice,  electron  collisionality  seems  to  be  much  smaller,  and  the  electron  velocity  dis¬ 
tribution  function  is  more  complex.  A  theory  for  the  partial  thermalization  case  is  presented 
in  Annex  B.  In  the  conclusions  section  of  that  Annex  we  comment  on  the  pending  work  and 
the  difficulties  of  implementing  that  theory  in  HPHall. 

3.4  Sheath  model  for  ideal  electron  emitter 

Active  electrodes  are  the  preferred  method  when  designing  two-stage  SPT-type  Hall  thrusters. 
Preferably  this  intermediate  electrode  is  made  of  a  high  electron  emissivity  material. 

This  new  element  inside  the  chamber  requires  proper  modelling  and  in  particular,  the 
plasma-wall  interaction  problem  is  modihed  with  respect  to  the  ceramic  wall  case  treated  in 
[10],  [2]  and  in  the  previous  section.  The  model  we  are  implementing  considers  an  active  elec¬ 
trode  operating  continuously  in  the  charge  saturated  regime(CSR)  due  to  its  high  electron 
emissivity.  The  potential  of  the  electrode  should  be  a  known  parameter,  whereas  the  electric 
current  through  the  electrode  would  be  the  output.  However,  at  present,  and  because  of  the 
unsteady  behavior  of  the  discharge  potential,  we  are  taken  the  electrode  current  as  the  known 
(and  constant)  parameter,  and  the  electrode  potential  as  the  unknown. 

Regarding  the  mathematical  formulation  of  the  model  we  can  expect  important  similarities 
with  respect  to  the  ceramic  wall  model,  including  the  usage  of  a  truncated  maxwellian  distri¬ 
bution  inside  the  sheath  for  primary  electrons.  In  fact,  the  governing  equations  are  the  same 
except  for  a  fundamental  difference:  the  zero  current  condition  is  no  longer  applicable. 

In  the  zero  Debye-limit  the  sheath  satishes  the  following  conservation  equations: 

-  ion  flux  conservation:  g^i  =  UiUri  =  const  =  g^q 

-  primary  electron  flux  conservation:  grp  =  UpUrp  =  ccmst  =  grpw 

-  secondary  electron  flux  conservation:  grs  =  UgUrs  =  const  =  grsw 

-  ion  energy  conservation:  (l/2)miU^j  -|-  e0  =  ccmst 

-  secondary  electron  energy  conservation:  {ll2)meulg  —  ecp  =  const  ~  —ecfcw 

-  primary  electron  Maxwell-Boltzmann  equilibrium: 

KT,  )  l  +  erfX^) 


np(0)  =  npQ  exp 
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-  primary  electron  semi-maxwellian  wall  flux; 
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where  erf{x)  is  the  error  function  and  subindexes  e,  p  and  s  refer  to  the  thermalized  electron 
population  in  the  bulk  of  the  plasma,  the  primary  electrons  and  the  emitted  electrons  respec¬ 
tively  (we  use  the  same  subindex  as  for  secondary  electron  due  to  similarity  of  the  formulation). 
Additionally  the  following  conditions  are  also  applicable  in  the  sheath  edge: 

-  conservation  of  electron  flux:  PreQ  =  QrpQ  —  drsQ 

-  emitted  charge  current:  jrQ  =  e{griQ  -  greo) 

-  plasma  quasi-neutrality:  riiq  =  UeQ  =  -|-  Usq 

-  charge  saturation  condition  (ideal  emitter);  U^ecpq)  =  Ui^ecpw) 

-  Bohm  condition;  f/(e0)"|Q  =  17(e0)"|^ 

where  U {ecj))  is  the  so-called  Sagdeev 's  potential,  a  function  that  appears  typically  in  sheath 
problems  as  a  solution  to  Poisson  A  equation.  Because  of  the  similarities  with  the  ceramic  wall 
sheath  the  form  of  this  function  is  the  same  as  in  the  problem  presented  in  the  previous  section. 

This  closed  set  of  equations  can  be  solved  in  a  non-dimensional  form  and  yields  results  in 
terms  of  two  main  parameters,  the  ion  Mach  number  and  the  non-dimensional  potential  sheath 
fall.  Figure  3.8  presents  the  dimensionless  sheath  magnitudes  for  the  active  electrode.  Both 
the  non-dimensional  current  to  the  electrode  sheath  and  the  relationship  between  the  emitted 
electron  flux  and  net  electron  flux  coming  from  the  sheath  are  shown,  for  sonic  and  supersonic 
ion  fluxes  (see  [2]).  The  higher  the  potential  fall  in  the  wall,  the  higher  is  the  electric  current. 
Notice  that  moderate  sheath  potential  falls  are  enough  to  emit  a  signihcant  electric  current. 


Figure  3.8:  Net  electron  wall  flux  to  electron  emitted  flux  ratio  (left)  and  non-dimensional  wall 
current  (right)  as  functions  of  the  non-dimensional  sheath  potential  fall  for  a  high-emission  electrode. 

In  chapter  6  this  model  is  used  to  simulate  an  intermediate  electrode  inside  a  two-stage  Hall 
thruster. 


Chapter  4 

IMPROVEMENTS  IN  PIC  SUBCODE 


This  chapter  explains  the  modihcations  implemented  in  the  PIC  sub-code  of  HPHall-2.  Even 
though  the  basis  of  the  method  are  mainly  conserved,  several  improvements  have  been  intro¬ 
duced  in  order  to  increase  the  accuracy  of  the  numerical  methods  and  thus  of  the  simulation. 
Some  of  these  changes  are  related  to  the  so-called  Bohm  condition  at  the  sheath  edges,  a  topic 
that  has  received  important  attention  lately  since  it  affects  greatly  overall  thruster  performance 
predictions.  Additionally,  other  topics  such  as  multiple  ionization  have  been  reviewed  and  new 
and  better  algorithms  have  been  implemented. 


4.1  Independent  doubly-charged  ion  population 

As  the  discharge  voltage  increases,  the  electron  temperature  increases  thus  augmenting  the 
multiple  ionization  fraction.  The  novelty  reported  here  is  to  treat  doubly-charged  ions  as  an 
independent  species  from  singly-charged  ions.  This  allows  to  have  reliable  weighting  magnitudes 
for  these  particles.  The  penalty  is  an  increase  of  a  20%  in  overall  computational  time  although, 
for  a  typical  case,  the  total  time  remains  within  the  order  of  an  hour  for  a  2.4GHz  PC. 

A  non-dimensional  parameter  that  measures  the  relative  importance  of  doubly-charged  ions 
is  the  local  ion  charge  number 

^i-\ — h 

which  usually  takes  values  between  1  and  1.1.  Figure  4.1  shows  the  effect  of  including 
double  ionization  for  a  SPT-lOO-like  HET  at  nominal  operation  conditions(300V).  Even  for 
this  moderate  voltage  the  anode  efficiency  r]a  increases  a  5%,  from  0.378  to  0.397,  mainly  due 
to  a  higher  propellant  utilization.  Notice  how  the  charge  number  increases  with  the  temperature 
from  the  anode  until  the  maximum  Tg.  Downstream  of  this  point,  the  charge  number  value  is 
due  to  ion  convection.  Additionally,  the  charge  utilization  is  rjq  =0.993.  Hofer  and  Gallimore 
obtained  experimentally,  for  the  NASA-173Mv2  thruster  and  300V,  a  lower  value,  rjq  =0.986, 
which  means  larger  multiple  ionization.  Thus,  the  need  to  consider  multiply  charged  ions  is 
justihed  even  at  moderate  voltages. 
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z(m) 


z(m) 


Figure  4.1;  Singly-  and  doubly-charged  ion  densities  along  the  channel,  and  nj++  respectively,  (left 
figure)  and  electron  temperature  Te  and  equivalent  ion  charge  number  Z*  along  the  channel  (right 
figure).  All  variables  are  evaluated  at  the  chamber  median,  r=  0.0425m. 


4.2  Cell  population  control 

4.2.1  Modified  ionization  algorithm 

PIC  methods  obtain  macroscopic  variables  (e.g.  densities,  particle  fluxes)  by  weighting  macro- 
paticle  magnitudes  to  the  nodes  of  the  computational  domain.  However,  in  order  to  avoid 
statistical  oscillations,  the  number  of  macroparticles  per  cell  must  be  high  enough.  A  good 
trade-off  between  accuracy  and  computational  cost  is  around  20  to  30  macroparticles  per  cell. 

HPHall-2  fails  to  maintain  an  acceptable  minimum  number  of  macroparticles  per  cell,  mainly 
near  the  anode  where  ionization  is  not  strong  enough  and  plasma  density  is  small.  This  fact 
has  important  consequences  since  thruster  efficiency  depends  on  energy  losses  to  walls  that  in 
turn  depend  on  ion  macroscopic  variables,  mainly  density  and  particle  flux  to  the  walls.  Thus, 
we  decided  to  modify  the  ionization  algorithm  in  order  to  control  better  the  ion  population. 
This  algorithm  is  still  of  the  MCC  type  and  ensures  a  minimum  number  of  particles. 

The  singly  charged  ion  mass  to  be  introduced  into  cell  (j,  k)  at  a  given  time  step  is  computed 
from  the  PIC  sub-code  (see  [3])  as: 


(4.2) 


where  is  the  cell  volume  and 


iPi)  jk 


iPi) jk  T  iPi)j+l,k  T  iPi)j,k+l  T  i^^i)j+l,k+l 
4 


(4.3) 


is  the  average  ionization  rate  for  cell  {j,  k). 

MCC  methods  are  based  on  performing  a  number  Nprobjk  of  probability  tests  to  determine 
the  real  number  of  macroions  to  be  created,  Nj^,  in  a  particular  cell  and  time  step.  The 
relationship  between  Nprobjk  and  {Ami)jk  is  given  by: 


Jk^probJkPjk  (^^^^i)jk') 


(4.4) 
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where  is  the  mass  of  the  created  macroions  and  Pjk  is  a  probability.  Once  Nprob,jk  is 

known  the  real  number  of  macroions  to  be  created  in  cell  (j,  k),  Njk,  is  determined  comparing 
Nprobjk  random  numbers  with  Pjk-  Notice  that  Nprobjk  represents  the  maximum  number  of 
macroions  that  could  be  created  in  the  cell  in  the  time  step  considered. 

Two  of  the  three  magnitudes  of  the  left-hand-side  of  equation  (4.4)  can  be  chosen  arbitrarily. 
HPHall  and  HPHall-2  take  Pjk  =  0.30—0.85  (depending  on  the  ionization  rate)  and  fix  either  the 
value  of  Mijk  or  Nprobjk  equal  to  one.  In  both  cases,  few  particles  are  generated  in  low  density 
regions.  In  order  to  control  better  the  number  of  particles  per  cell,  the  ionization  algorithm 
has  been  improved  in  the  following  way.  First,  a  law  for  the  probability  Pjk  is  chosen;  it  can  be 
constant  (e.g.,  0.5)  or  a  function  of  the  ionization  rate  as  in  the  existing  codes.  Second,  lower 
and  upper  limits  of  macroions  per  cell  are  dehned,  N^in  and  N^ax-  Then,  three  situations  are 
considered,  depending  on  the  numbers  of  existing  macroions  in  the  cell  Nijk'- 

1.  If  Nijk  <  Nfnin  then  NprobjkPjk  =  Nmin  —  ^ijk  (what  eusures  that  the  minimum  of 
macroions  is  respected)  and  Mijk  is  obtained  from  equation  (4.4). 

2.  If  Nmin  <  P^i,jk  <  ^max  then  Mj jfc  is  choseu  close  to  the  average  value  of  the  mass  of  the 
macroions  in  the  cell  and  the  integer  Nprobjk  is  computed  from  equation  (4.4). 

3.  If  Nijk  >  Nmax  then  Npj-objk  =  1  (what  reduces  the  computational  cost  and  ensures  that 
there  are  no  cells  with  too  many  particles)  and  Mijk  is  obtained  from  equation  (4.4). 

4.2.2  Results 

The  main  advantages  of  the  new  algorithm  with  respect  to  the  previous  one  are: 

1.  it  reduces  statistical  oscillations  by  ensuring  a  minimun  number  of  particles  in  every  cell. 

2.  the  computational  cost  is  reduced  by  controlling  the  maximum  number  of  macroions  per 

cell. 

3.  macroscopic  variables  are  weighted  properly  near  the  walls  since  an  acceptable  number 
of  macroions  exists. 

Figure  (4.2)  shows  the  differences  between  HPHall  algorithms  results  and  the  results  ob¬ 
tained  with  the  algorithms  proposed  here.  Clearly,  the  aforementioned  advantages  are  achieved. 
The  new  algorithm  is  capable  of  limiting  the  ion  population  per  cell  between  20  and  50,  ap¬ 
proximately.  In  comparison,  HPHall  generated  too  many  particles  at  the  exit  of  the  chan¬ 
nel  (~  130),  what  increased  the  computational  cost  without  improving  simulation  accuracy, 
whereas  it  worked  with  few  particles  (<  10)  near  the  anode,  what  caused  large  temporal  oscil¬ 
lations,  in  the  scale  At,  of  statistical  character.  Indeed,  HPHall  added  an  artihcial  background 
density  of  10^^  (part/m^)  to  Ue  to  avoid  ’near-singular’  values.  These  minima  are  clearly  observ¬ 
able  in  the  plots  of  ne{t)  for  HPHall.  The  new  population  control  algorithm  solves  excellently 
that  issue.  The  statistical  oscillations  of  (j){t)  follow  those  of  ne{t).  Differences  on  the  temporal 
evolution,  in  scales  much  larger  than  At  are  due  to  other  improvements  of  the  code  [13,  14]. 
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Figure  4.2;  Comparison  between  HPHall-2  (left  column)  and  the  new  hybrid  simulation  code 
(right  column);  Contour  Number  of  macroions  per  cell,  temporal  evolution  of  plasma  density 
and  potential  in  a  node  of  the  anode  sheath  edge. 
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4.3  Bohm  condition  fulfilment:  surface  weighting 

References  [12]  and  [15]  showed  that  a  correct  weighting  of  plasma  magnitudes  at  the  boundaries 
of  the  quasineutral  computational  domain  was  a  delicate  issue.  In  this  section  we  review  the 
methods  that  have  been  used  in  HPHall  and  HPHall-2  and  we  present  a  new  scheme  which 
compares  advantageously  with  previous  ones.  The  influence  of  the  different  schemes  over  the 
temporal  oscillations  is  briefly  discussed  too. 

Figure  4.3  shows  the  thruster  geometry,  the  two  test  meshes  and  the  magnetic  held  prohle 
used  in  this  section  work.  Parameters  correspond  to  a  SPT-70  thruster.  To  facilitate  the 
understanding  of  wall  effects,  we  use  a  rather  simple  geometry  and  a  B-£eld  satisfying  =  0 
and  Bj.  (X  1/r;  in  this  way  V  ■  B  =  0,  but  V  x  B  7^  0. 


Figure  4.3;  Geometry,  test  meshes  and  B-field  for  this  work.  Typical  execution  times  (for  a  Pentium 
III  at  2.4  GHz)  to  reach  a  steady-state  response  are  included. 


4.3.1  Weighting  schemes 

Volumetric  Weighting(VW) 


This  is  the  original  scheme  used  by  Fife  [16].  With  this  scheme,  the  magnitudes  at  nodes  placed 
along  the  domain  boundaries  are  computed  using  the  typical  PIC  weighting  algorithm.  The 
ion  density  at  a  generic  node  T  is  determined  from 


(neT)vW 


ZpTnp 

St{Tpi  Zp), 
rrii 


(4.5) 


where  AVr  is  the  volume  associated  to  node  T,  p  refers  to  a  particle  being  within  the  volume 
of  inhuence  of  the  node,  rup  and  {zp,rp)  are  particle  mass  and  position,  Zp  is  one  for  single 
ions  and  two  for  double  ions,  and  ST{rp,  Zp)  is  a  bi-linear  type  weighting  function  for  node  T. 
This  algorithm  was  applied  to  boundary  nodes,  even  though  the  volume  of  inhuence  at  the 
boundaries  is  diherent  (see  Figure  6  of  Ref.  [15]).  The  one-side,  asymmetric  weighting  at 
the  boundary  nodes  tends  to  underestimate  [overestimate]  magnitudes  that  increase  [decrease] 
toward  the  boundary.  The  error  can  be  large  for  magnitudes  that  present  large  gradients  near 
the  sheath  boundary  (which  include  Ue,  Vj,  and  0,  but  not  jj. 

In  [12]  we  demonstrated  that,  using  VW  and  as  the  cell  size  of  the  mesh  is  reduced,  the  ion 
how  at  the  boundary  increases  significantly  and  larger  plasma  gradients  over  the  whole  radial 
section  develop.  However,  for  the  mesh  sizes  the  ion  how  remains  well  subsonic  (i.e.  still  far 
from  satisfying  the  Bohm  condition),  making  it  evident  that  a  diherent  numerical  approach 
should  be  used  in  order  to  weight  correctly  magnitudes  at  the  boundaries. 
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Corrected  Volumetric  Weighting(CW) 

This  scheme  was  proposed  in  Ref.  [15].  The  idea  is  to  use  the  information  of  two  nodes  [the  one 
in  the  boundary(T)  and  the  internal  one  next  to  it(T+l)]  in  order  to  correct  the  0(Ar)  error 
induced  by  the  asymmetric  VW.  With  the  information  of  two  nodes,  the  error  can  be  reduced 
to  O(Ar^).  The  corrected  plasma  density  is, 

3  1 

(^er)cw  =  2(wr)vw  —  ^(ne^T+i)vw,  (4.6) 

with  (ne)vw  the  plasma  density  obtained  from  the  VW.  (ner)cw  is  smaller  than  (ner)vw  for 
usual  conditions.  This  correction  lies  on  internal  (i.e.  PIC-computed)  magnitudes  only,  but,  at 
the  same  time,  it  depends  on  a  second  node,  and  thus  errors  depend  on  the  cell  size.  Reference 
[15]  proved  the  CW  algorithm  to  give  a  reasonable  fulhllment  of  the  Rohm  condition  if  the 
mesh  is  hne  enough. 

Bohm  Condition  Forcing  (BCF) 

This  algorithm  was  presented  in  Ref.  [15].  It  consisted  on  forcing  the  ion  flux  to  satisfy  the 
(simple  form  of  the)  Bohm  condition 

j*  ■  N  >  ne\/ZTe/mi  (4.7) 

whenever  this  was  not  achieved  naturally  by  the  VW  algorithm;  Z  is  the  average  ion  charge- 
number.  This  the  BCF  algorithm  for  the  density  is 

{neT)BCF  =  min  |  (neT)vw  1  (4.8) 

[e^ZTe/nii  J 

With  this  algorithm,  UeT  is  lower  in  most  points  and  instants  than  in  the  VW  and,  therefore, 
the  electric  held  is  larger  near  the  boundaries.  In  Ref.  [15]  we  showed  that  this  algorithm 
provides  excellent  results  even  for  relatively  coarse  meshes.  This  algorithm,  contrary  to  the 
previous  ones,  relies  on  a  property  dictated  by  the  space-charge  sheath  and,  hence,  extrinsic  to 
the  PIC  code. 


Surface  Weighting  (SW) 

Although  the  CW  algorithm  was  promising,  the  need  of  a  fine  mesh  was  evident.  The  idea  of 
SW  is  to  count  the  particles  that  cross  the  surface  of  influence  of  the  node  at  the  boundary. 
Thus,  SW  is  a  PIC-intrinsic  algorithm  but,  at  the  same  time,  is  independent  of  the  cell  size. 
The  density  Ue  in  the  boundary  node  T,  as  computed  by  SW,  is 


(^^er)sw 


1 

AStAI 


E 


ZpTTlp  1 

m  Vp-N’ 


(4.9) 


where  the  sum  extends  to  all  the  ions  that  have  crossed  the  boundary,  N  is  the  normal  to  the 
wall,  ASt  is  the  surface  of  influence  of  node  T,  and  At  is  the  timestep  used  by  the  PIC  code. 


4.3.  BOHM  CONDITION  FULFILMENT:  SURFACE  WEIGHTING 


37 


This  computation  for  UeT  can  be  used  whenever  the  flow  of  particles  is  only  in  one  direction, 
like  in  the  boundary,  where  no  ions  are  coming  back.  The  obtained  (neT)sw  is  a  time  average  in 
At  of  HeT,  and  its  error  is  only  O(Ar^)  (associated  to  the  discretization  of  the  electric  potential, 
mainly) . 

Similarly,  the  ion  current  density  at  the  boundaries  is  given  by 


(jjr)sw  - 


e 

ASrAt 


E 


ZpTflp  Vp 

mi  Vp-N’ 


(4.10) 


4.3.2  Comparison  of  profiles  for  the  different  schemes 

Figure  4.4  shows  some  radial  prohles  obtained  with  different  algorithms  and  mesh  sizes.  The 
new  SW  is,  like  CW,  intrinsic  to  the  PIC  model,  but  at  the  same  time  provides  reasonably  good 
results  even  for  coarse  meshes.  The  presheath  potential  drop  for  SW  and  mesh  1  is  around 
75%  of  the  potential  drop  achieved  with  BCF,  and  the  difference  decreases  to  5%  just  by 
dividing  by  two  the  mesh  size  around  the  walls.  This  large  difference  is  mainly  due  to  the  steep 
slope  that  exists  in  the  density  and  potential  prohles  near  the  boundaries  of  the  quasineutral 
domain.  Other  magnitudes,  like  the  ion  density  current  how  to  the  walls  (j^-p)  does  not  change 
as  dramatically  (see  hgure  4.5).  The  radial  potential  drop  is  not  adequate  to  estimate  the 
numerical  errors  induced  by  the  hybrid  model  because  it  is  very  sensitive  to  the  values  near 
the  boundaries,  where  the  slope  is  extremely  high.  Magnitudes  with  hatter  slopes  are  probably 
more  suitable  to  analyze  the  accuracy  of  the  schemes  (as  it  is  the  case  of  jj'p).  However,  radial 
prohles  and  their  development  give  a  good  understanding  of  the  problem,  even  if  the  exact 
values  at  the  boundary  are  too  sensitive  to  take  them  into  account.  Using  as  a  criterium, 
we  probably  can  say  that  BCF  and  SW  are  more  accurate,  since  the  change  in  ion  how  with 
the  mesh  size  is  not  as  high  as  it  is  in  CW. 


Figure  4.4;  Radial  profiles  of  electrostatic  potential  in  z  =  0.0256  m.  (a)  Comparison  of  results  with 
the  different  algorithms  in  mesh  1.  (b)  Comparison  of  results  with  the  different  algorithms  in  mesh  2. 

As  observed  in  Ref.  [12],  the  ion  how  to  the  walls  increases  when  the  cell  size  decreases. 
Also,  the  schemes  that  provide  more  accuracy  (BCF  and  SW)  give  a  bigger  ion  how  to  the 
walls.  The  ion  how  can  even  be  double  with  BCF  and  SW  than  it  is  with  CW.  Although  BCF 
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z  (m)  z  (m)  z  (m) 

Figure  4.5:  Axial  profiles  of  ion  current  density  to  the  inner  wall  for  (a)  mesh  1  and  (b)  mesh  2.  In 
figures  (a)  and  (b),  the  dashed  line  is  the  theoretical  Bohm  flux  (jBohm)  for  SW  and  the  corresponding 
mesh,  (c)  Comparison  of  sensitivities  to  the  mesh  size  of  CW,  SW  and  BCF.  Dashed  lines  represent 
results  for  mesh  1  and  solid  lines  represent  results  for  mesh  2. 

and  SW  have  similar  accuracy  -  if  we  consider  the  sensitivity  to  mesh  size  a  measure  of  the 
error  there  is  a  difference  in  the  ion  flow  at  the  exit  of  the  channel.  This  difference  is  small 
compared  to  the  difference  between  those  schemes  and  CW  or  VW,  but  it  is  still  important 
and  remains  to  be  studied. 

In  hgures  4.5(a)  and  4.5(b)  the  resulting  ion  flow  with  SW  is  compared  to  the  simple  Bohm 
condition  (4.8).  With  mesh  1,  both  results  differ  greatly,  probably  due  to  the  very  different 
values  of  Ue  at  the  boundaries  (we  have  already  commented  the  sensitivity  of  boundary  values). 
The  comparison  is  much  better  for  mesh  2,  where  both  values  are  practically  the  same. 

To  sum  up,  surface  Weighting  and  Bohm  Condition  Forcing  yield  the  best  results.  The 
Corrected  Weighting  requires  too  hne  of  a  mesh,  although  it  can  be  good  to  conhrm  the  results 
given  by  the  two  other  methods.  Between  SW  and  BCF,  SW  has  the  advantage  of  being  more 
self-consistent  since  it  is  PIC-intrinsic.  If  the  error  is  estimated  from  the  sensitivity  of  the 
current  to  the  wall  with  respect  to  the  mesh  size,  BCF  and  SW  have  similar  performance.  SW 
hts  the  Bohm  condition  (4.8)  for  reasonable  mesh  sizes,  which  conhrms  its  accuracy. 

4.3.3  Discharge  oscillations 

One  of  the  uses  of  simulation  codes  is  the  study  of  plasma  oscillations  in  Hall  thrusters.  There¬ 
fore,  the  effect  of  the  different  schemes  on  time  oscillations  is  an  important  criterium  to  decide 
how  appropriate  the  different  schemes  are.  In  hgure  4.6,  we  present  some  time  prohles  for  the 
discharge  current.  Two  well-defined  frequencies  appear;  one  around  15  kHz,  which  seems  to  be 
related  to  the  breathing  mode,  and  another  one  near  170  kHz,  which  is,  approximately,  the  ion 
transit  frequency.  These  frequencies  do  not  change  with  the  mesh  size  or  the  weighting  scheme, 
but  their  amplitude  depends  greatly  on  these  parameters.  The  15-kHz-oscillation  amplitude 
is  highly  dependent  on  the  weighting  scheme:  for  mesh  2,  the  amplitude  jumps  from  35%  for 
SW  to  15%  for  BCF.  It  happens  similarly  for  mesh  1.  In  the  case  of  BCF,  there  must  be  two 
modes  very  near  15  kHz  that  produce,  when  they  are  added,  the  variable  amplitude  seen  in 
hgures  4.6(b)  and  4.6(d)  (it  is  especially  clear  in  this  last  one).  On  the  other  hand,  the  170  kHz 
oscillation  is  not  as  dependent  on  the  weighting  scheme.  For  mesh  2,  its  amplitude  is  15%  of 
the  average  current  with  both  BCF  and  SW.  For  mesh  1,  the  amplitude  decreases  to  8%  when 
we  are  using  SW  and  is,  again,  around  15%  with  BCF.  BCF  seems  to  be  able  to  reproduce 
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this  mode  correctly  for  coarse  meshes,  whereas  with  SW  the  amplitude  of  the  high  frequency 
oscillations  is  rather  dependent  on  mesh  size. 

To  sum  up,  BCF  perturbs  less  the  170  kHz  mode,  but  it  provokes  a  variable  amplitude  for  the 
breathing  mode,  perhaps  due  to  the  superposition  of  two  modes  of  very  similar  frequency.  SW 
is  not  as  effective  reproducing  the  high  frequency  oscillations  as  BCF  is,  but  with  hue  enough 
mesh  SW  yields  the  same  amplitude.  There  are  still  important  differences  in  the  breathing 
mode  between  both  schemes,  which  remain  to  be  studied. 


t  (|XS)  t  (|XS) 


Figure  4.6;  Time  evolution  of  discharge  current  for  different  schemes  and  mesh  sizes. 
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4.4  PIC  subcode  accuracy  improvement 

This  work  intends  to  correct  certain  inconsistencies  found  in  both  HPHall  and  its  upgraded 
version  HPHall-2  that  mainly  affect  the  PIC  subcode.  The  problems  are  related  to  the  ‘leapfrog’ 
integration  scheme  [17]  and  its  connection  with  the  weighting  algorithms  employed  to  transform 
macroparticle  information  into  macroscopic  variables  such  as  densities,  particle  fluxes,  energy 
fluxes,  etcetera. 

Macroparticle  equations  of  motion  are  deduced  from  Newton 's  law; 

=  F{xp,Vp,t)  =  qp(^E{xp,t)  +  VpX  B{xp)^  (4.11) 

where  Mp,  Xp  and  Vp  are  respectively  macroparticle  charge,  mass,  position  and  velocity  and 
E  and  B  are  the  electric  and  the  magnetic  helds  respectively.  Besides,  another  force  term  can 
be  added  accounting  for  both  ion-neutral  and  charge-exchange  collisions.  The  ‘leapfrog’  scheme 
permits  to  integrate  numerically  the  equations  of  motion  of  all  the  macroparticles  using  the 
following  algorithm; 


^p,t+At/2  At/2  T  '  F(^Xp  ^R^  (4.12) 

^p,t+At  ^p,t  T  '  '^p,t+At/2 

where  At  is  the  PIC  integration  time  step,  F  is  the  force  on  the  macroparticle,  and 
subindexes  stand  for  the  time  step  in  which  variables  are  computed.  Notice  that  position 
and  velocity  are  computed  at  instants  that  differ  At/2,  hgure  4.7. 


At  in  X  At  in  X 


t  -  At/2  t  +  At/2 


At  in  V 


Figure  4.7;  ‘Leapfrog’  integration  scheme  for  macroparticle  motion. 


The  numerical  error  associated  to  the  ‘leapfrog’  scheme  is  O(At^)  in  velocity  and  O(At^) 
in  position,  provided  that  forces  acting  on  macroparticles  are  known  with  O(At^)  accuracy. 
Therefore,  the  electric  held,  which  is  computed  by  the  electron  subcode,  must  have  O(At^) 
accuracy.  In  turn  the  electron  sub-code  needs,  from  the  PIC  sub-code,  the  electron  density, 
TTg,  and  the  ion  current  hux,  j/.  Thus,  for  the  electric  held  to  be  known  with  0(A/^)  accuracy, 
Ue  and  ji  must  be  computed  with  the  same  precision.  Both  HPHall  and  HPHall-2  fail  to  keep 
this  level  of  accuracy.  We  explain  next  the  corrections  made  to  the  diherent  macroparticle 
algorithms  in  order  to  solve  this  precision  problem. 
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4.4.1  Time-centered  weighting  of  macroscopic  variables 

The  ion  density,  Ui,  and  the  ion  current  flux,  ji,  at  a  generic  node  (j,  k),  are  obtained  by  the 
PIC  sub-code  from  the  weighting  formulas[3] 

{Pi)jk  =  A-tr  ^  MpSjk{Xp),  (4-14) 

7/7'^^  V  weight,  ^ 

{ji)jk  =  ^  '^MpVpSjk{xp),  (4-15) 

fTT'i^^weightjjk  ^ 

where  Sjk{x)  is  the  (bi-linear)  weighting  function,  AVweight,ifc  volume  of  influence  of  the 

node[18],  the  summation  extends  to  all  macroions  in  the  vicinity  of  the  node  (where  Sjk  ^  0), 
Mp  is  the  macroparticle  mass  (which  varies  as  a  result  of  ionization),  and  e  and  are  the 
charge  and  mass  of  the  ion  species.  Similar  weighting  formulas  apply  for  the  pressure  and 
higher-order  variables. 

In  HPHall-2  the  ionization  process  is  modelled  differently  for  the  macroneutrals  than  for  the 
much  lighter  (singly-charged)  macroions;  in  each  timestep,  ionization  leads  to  the  loading  of 
new  macroions  in  the  domain,  but  to  a  weak  reduction  of  the  mass  of  each  macroneutral  [16]. 
A  similar  process  regulates  the  creation  of  doubly-charged  macroions  from  macroneutrals  or 
singly-charged  macroions.  Here,  for  simplicity,  we  consider  that  only  one  species  of  (singly- 
charged)  ions  is  present,  so  that  mass  variation  affects  to  macroneutrals  only. 

In  HPHall-2,  the  macroscopic  magnitudes  are  computed  using  position  and  mass  at  time  t 
and  velocity  at  time  t  —  At/2,  so  their  precision  is  only  0{At).  In  order  to  obtain  all  weighted 
variables  at  instant  t  —  At/ 2  with  precision  O(At^),  the  position  and  mass  of  macroparticles  at 
time  t  —  At/ 2  must  be  used.  The  simple  interpolation  rules 


t— Ai/2 

(4.16) 

2 

AIp^t-At/2  = 

^p,t  + 

(4.17) 

2 

provide  precision  O(Af^).  Besides,  the  velocity  obtained  with  the  ‘leapfrog’  scheme  has  a 
numerical  error  of  the  same  order.  Therefore,  the  weighting  process  has  the  required  0{AU) 
accuracy. 

4.4.2  Time-centered  loading  of  particles 

At  each  timestep,  HPHall-2  loads  macroparticles  into  the  simulation  domain  in  order  to  account 
for  different  processes,  such  as  neutral  gas  injection  and  ionization,  particle  wall  recombination, 
etcetera.  Care  must  be  taken  regarding  the  instant  of  loading.  The  MCC  method  of  HPHall-2 
provides  an  statistical  distribution  of  the  location  and  velocity  of  the  loaded  particles,  but  not  on 
the  instant  of  loading.  Thus,  the  location  and  velocity  are  implicitly  assigned  to  instants  t  and 
t  —  At/2,  respectively.  The  enhancement  proposed  here  is  to  determine  randomly  the  instant  of 
loading,  that  is  f  —  fAt,  with  /  G  [0, 1].  Thus,  the  location  and  velocity  provided  by  the  MCC 
method  correspond  to  that  instant.  It  is  then  necessary  to  center  the  position  and  velocity  of 
the  loaded  particles  at  times  t  and  t  —  At/2,  with  errors  ©(At^^)  and  O(At^),  respectively,  in 


42 


CHAPTER  4.  IMPROVEMENTS  IN  PIC  SUBCODE 


order  to  be  consistent  with  the  leapfrog  scheme.  The  algorithm  of  Cartwright  et  a/.  [19]  was 
adapted  to  HP  Hall  in  Ref.  [20].  However,  this  complex  procedure  is  not  completely  necessary 
since  loaded  particles  are  always  macroneutrals  (injection,  wall  accomodation  or  chamber  flux 
modelling)  and  it  is  possible  to  assume  constant  mass  and  constant  speed  for  the  time  step 
in  which  the  particle  is  loaded.  However,  in  order  to  be  consistent  with  the  leapfrog  scheme 
position  and  velocity  must  be  computed  at  t  and  t  —  At/2  respectively  from  position  and  velocity 
at  the  instant  of  injection. 

4.4.3  Comments  and  results 

The  previous  modihcations  have  required  several  changes  in  the  code  since  now  it  is  necessary 
to  distinguish  two  instants,  t  —  At/2  and  t.  A  second  detail  to  take  into  account  is  that  particles 
that  are  outside  the  computational  domain  at  t  cannot  be  removed  from  the  computations  until 
weighted  magnitudes  are  computed  at  t  —  At/2.  And  inversely,  injected  particles  can  be  inside 
the  domain  at  t  but  outside  at  t  —  At/2  and  thus,  those  particles  must  not  be  accounted  for 
when  weighting  variables  at  t  —  At/2.  This  is  particularly  important  in  order  to  reproduce  the 
desired  mass  flow  through  the  injector. 

However,  there  is  a  third  improvement  that  is  lacking  in  the  current  implementation.  It  is 
related  to  the  computation  of  the  exact  instant  in  which  a  particle  impacts  the  wall  and  the 
associated  variables  (mass,  position  and  velocity). 

The  main  consequence  of  these  improvements  is  an  increase  of  the  PIC-subcode  accuracy. 
This  increase  in  time-accuracy  of  the  code  has  some  relevant  effects  on  the  simulation  results 
regarding  the  mass  conservation.  The  mass  flow  along  the  channel  is  conserved  more  accurately. 
Figure  4.8  depicts  ion,  neutral  and  total  mass  flow  for  an  SPT-100  like  simulated  thruster  (300 
V,  4.8  mg/s)  computed  using  both  HPHall-2  and  the  modihed  algorithms.  HPHall-2  presents 
a  mass  flow  loss  of  1.5%  between  the  anode  and  the  channel  exit,  while  the  error  is  below  0.4% 
with  the  new  method.  Note  that  the  most  important  change  takes  place  in  the  ion  flow.  This 
is  likely  due  to  proper  account  for  the  macroparticle  position  at  (t  —  At/2). 


Figure  4.8:  Mass  flows  in  a  SPT-100  like  simulated  Hall  thruster  computed  by  HPHall-2  (dashed 
lines)  and  the  new  algorithms  (solid  lines). 


Notice  that  both  results  have  been  obtained  running  the  program  for  more  than  20000  itera¬ 
tions  and  therefore  non-stationary  effects  are  removed  out  when  averaging  discharge  properties. 


Chapter  5 

NEW  ELECTRON  SUBCODE 


This  chapter  is  devoted  to  the  electron  subcode.  The  new  electron  differential  formulation  is 
explained  in  depth  and  the  main  advantages  are  highlighted.  An  additional  section  for  the 
description  of  the  near  anode  region  is  also  included  in  order  to  explain  the  closure  of  the 
electron  equations  in  this  region  of  the  thruster. 


5.1  One-dimensional  differential  formulation 


5.1.1  Magnetic- field-based  reference  frame 

Electrons  are  highly  magnetized  inside  Hall  thrusters  due  to  the  strong  magnetic  held.  Thus 
electron  motions  parallel  and  perpendicular  to  the  magnetic  held  are  completely  diherent.  This 
fact  makes  it  convenient  to  describe  electron  physics  in  a  magnetic-held-based  reference  frame. 
The  induced  magnetic  held  is  negligible  with  respect  to  the  externally  applied  magnetic  held. 
Furthermore,  the  latter  is  assumed  irrotational  (V  x  R  =  0)  and  solenoidal  (V  ■  R  =  0) 
and  therefore,  potential  {a{r,z))  and  stream  (A(r,  2:))  functions  exist  (see  hgure  1.1).  These 
functions  are  dehned  by  the  following  systems  of  equations; 


dr 


Br 


dz 


R. 


A  : 


d\ 

dr 


—rB 


Zl 


(5.1) 


The  resulting  reference  frame  is  dehned  by  the  coordinate  set  {a,  6,  A)  and  the  corresponding 
set  of  unit  vectors  {u\\,U0,  u±)  where  ug  is  a  unit  vector  along  the  azimuthal  direction,  -un  =  B/B, 
u±  =  u\\  X  Ug,  B  is  the  magnetic  held  vector  and  R  is  its  modulus.  The  diherent  vectors  and 
coordinates  are  shown  schematically  in  hgure  5.1. 

Additionally,  it  is  useful  to  dehne  spatial  coordinates  for  arc  lengths  along  the  streamlines 
(A  =  const)  and  along  the  lateral  walls  (r  =  r{z))  called  y  and  s  respectively: 


da  X  R  ’ 


ds  1 

d\  Hz)  r{B-N)' 


(5.2) 


with  N  a  vector  perpendicular  to  the  walls  and  oriented  towards  increasing  radii,  hg.  5.1  . 
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Anode 


Magnetic  -  >  ^ 

streamline  \  '^r. 

u,;- 


N 


n» 


Outer  wall 


Inner  wall 


Figure  5.1;  Magnetic-field-based  reference  frame. 


5.1.2  Electron  equations.  General  hypothesis. 

Electron  governing  equations  are  charge,  momentum  and  energy  conservation.  In  the  following 
subsections  each  one  of  them  will  be  analyzed  separately. 

Charge  conservation  equation 

Charge  conservation  law  can  be  written  as 


V-J=0  (5.3) 

where  j  =  ji  +  is  the  total  charge  current  composed  of  the  electron  current,  je  =  —eUeUe, 
and  the  total  ion  current  ji  =  enj+'Uj+  -|-  2enj++'Uj++.  Notice  that  Uj  stands  for  the  j-th  species 
velocity  vector. 

Momentum  conservation  equation 

Electron  motion  can  be  approximated  as  diffusive  (luel  ^  ^/Te/mg)  because  of  the  electromag¬ 
netic  forces.  Thus,  inertial  terms  can  be  neglected  in  the  momentum  conservation  law  which 
yields  the  following  simplihed  equation; 

0  =  — VUeTe  —  ene{E  +  Ue  B)  —  Uemel^eUe,  (5.4) 

where  rUe  and  Tg  are  the  electron  mass  and  temperature,  respectively,  E  is  the  electric  held 
vector  {E  =  —  V0),  cf)  is  the  electric  potential  and  Pg  =  Pen  +  ^wm  +  ^ano  is  the  total  electron 
collision  frequency.  Notice  that  electron-neutral  collisions  {Pen)  dominate  with  respect  to  other 
types  of  binary  collisions[21],  but  two  extra  terms  are  added  in  order  to  account  for  near- wall 
conductivity  (p^m)  and  anomalous  diffusion  {Pano)-  A  detailed  model  of  these  terms  will  be 
given  later. 

The  collision  term  is  negligible  with  respect  to  the  other  terms  along  the  streamlines  and, 
thus,  projecting  equation  (5.4)  onto  (un)  yields; 


0 


d  .  rn\  94> 

-{u,T,)+en,- 


(5.5) 
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The  temperature  is  approximately  constant  along  the  streamlines  (as  we  will  justify  later), 
and  therefore,  equation  (5.5)  can  be  integrated  to  yield  the  Maxwell-Boltzmann  equilibrium 
relationship; 

e0(a.  A,  t)  =  Te(A,  t)  In  +  e0*(A,  t)  (5.6) 

where  0*  is  the  so-called  thermalized  potential,  whose  dehnition  depends  on  the  chosen  reference 
density  n*,  that  can  be  dehned  as  a  constant  or  as  an  average  density  along  each  streamline. 
The  projection  of  equation  (5.4)  onto  the  azimutal  and  perpendicular(-ux)  directions  yields: 


0 


0  =  -eUeU^eB  -  nemeVeUBe  , 


(ueTe)  -f-  eUerB 


(90 


4-  eUeUeeB  -  nemeVeU^e  ■ 


(5.7) 

(5.8) 


The  Hall  parameter  is  dehned  as  the  ratio  between  the  electron  gyro-frequency  and  the  electron 
collisional  frequency,  0  =  cUe/zze  =  {^B) / {meVe)-,  and  its  value  is  higher  than  100  for  typical  Hall 
thruster  conditions  (0  3>  1).  If  expressions  (5.7)  and  (5.8)  are  combined  with  the  Maxwell- 
Boltzmann  equilibrium  relationship  and  the  Hall  parameter  dehnition,  the  following  expressions 
can  be  derived  for  the  electron  velocities  u±e  and  uoe'. 


UOe  fAu_i_e  1 


r(3  /  (90 
1  +  02 


eue  d\ 


(5.9) 


Equations  (5.9)  show  respectively  the  closed-drift  described  by  electrons  and  a  generalized 
Ohm’s  law  relating  the  electron  velocity  with  the  electric  held  and  the  pressure  gradient. 


Energy  conservation  equation 

Energy  conservation  can  be  expressed  in  the  dihusive  limit  as: 

+  V  ■  QueTeUe  +  =  je  '  E  -  heQ(i(Te)Ei  ,  (5.10) 

where  is  the  heat  conduction  hux,  ■  E  is  the  Joule  heating  and  neOii{Tf.)Ei  are  energy  losses 
due  to  ionization.  Here,  fie  is  the  ionization  rate,  ai{Te)  is  the  non-dimensional  net  energy  loss 
per  electron  and  Ei  is  the  Xenon  hrst  ionization  energy.  Notice  that  elastic  collision  energy 
losses  are  negligible  with  respect  to  ionization  losses[16]. 

The  ionization  rate  (he)  can  be  computed  as  a  combination  of  the  diherent  ionization  pro¬ 
cesses: 


<9  /3  _ 

—  I  -UeR 


dt\2 


He  =  n. 


'i-\- 


+  h, 


^e+ 
2+  + 


+  2h. 


i++ 


=  rheUniC^ 


+  2C. 


+  neUi+C,^ 


(5.11) 


where  is  the  neutral  density  and  variables  C,j  depend  on  the  electron  temperature [16].  These 
functions  are  plotted  in  hgure  5.2(a)  for  each  one  of  the  diherent  processes.  Notice  that  multi¬ 
ple  ionization  is  not  negligible  for  moderate-to-high  temperatures  what  makes  it  necessary  to 
account  for  doubly-charged  ions  in  Hall  thrusters.  Dugan’s  model  for  the  ionization  cost  factor 
Q;j(Te),  hg.  5.2(b),  is  used. 
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Figure  5.2: 
processes. 


Te/E,. 


(a)  Comparison  of  ionization  rate  parameters  Cj(^e)  for  Xenon  different  ionization 

(b)  Normalized  electron  production  cost  ai  versus  normalized  temperature  6  = 


The  appropriate  diffusive  model  for  heat  conduction  according  to  reference  [22]  is: 

(5/2)neTeVTe  +  eqeX  B  +  rUel'eqe  =  0 

Projecting  it  onto  the  magnetic  reference  frame  yields: 

5  UeTeP  dTe  5  UgrTe  13  dTe 

^""^“2  e  da  ’  “2  e  1  +  (3^  dX  ' 

The  estimate  of  the  relative  order  of  magnitude  of  the  derivatives  of  Tg  yields 

dTcldX  qfe  q\\c 


(5.12) 


Therefore,  for  any  plausible  ratio  between  heat  conduction  components,  the  variation  of  the 
electron  temperature  along  the  streamlines  seems  negligible,  justifying  the  approximation  of 
isothermal  streamlines,  i.e.  Te  =  Te(A,t). 


5.1.3  Evolution  equations  across  streamlines 

Let  us  consider  a  general  conservation  law  of  the  form,  V  ■  u  =  /,  for  certain  functions  u  and 
/.  The  integration  of  this  expression  between  two  different  streamsurfaces  located  at  A  and 
X  +  dX,  being  dX  an  inhnitesinial  quantity,  yields 


(5.13) 


where  the  summation  is  extended  to  the  inner  [i  =  1)  and  the  outer  {i  =  2)  walls,  n  points  out 
of  the  wall,  fig.  5.1,  and  r(A)  stands  for  the  streamline  located  at  A. 
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The  electron  subcode  of  HPHall-2  was  designed  for  ceramic  walls  only  and  thus  considered 
that  the  discharge  current  was  constant  along  the  channel.  This  simplihed  enormously  the 
electron  continuity  equation.  The  present  code  eliminates  that  constraint  by  accepting  net 
exchanges  of  electric  current  at  the  lateral  walls.  The  previous  integration  procedure  is  applied 
to  the  charge  equation  (5.3)  in  order  to  obtain  a  general  charge  conservation  equation  for  the 
electron  current  across  the  streamlines; 


j-n 


d\ 


d\ 


1=1,2 


B-N 


Ii,e{\)=  I  {ji,e-U±)  dS. 

ds{X) 


(5.14) 


Here  j  ■  n  is  the  electric  current  density  to  the  lateral  (sheaths  and)  walls,  and  the  ion  cur¬ 
rent  density  ji  is  computed  by  the  particle  sub-code.  This  charge  conservation  equation  is  one 
of  the  main  novelties  of  the  present  model.  SPT-type  thrusters  normally  have  ceramic  walls 
and  therefore  no  net  charge  flux  to  lateral  walls  exists  {j-n  =  0).  However,  some  new  de¬ 
signs  incorporate  intermediate  segmented  electrodes[23]  in  the  lateral  walls  in  order  to  enhance 
performances  and  control.  This  new  formulation  permits  to  simulate  these  conhgurations.  Ad¬ 
ditionally,  a  more  physically  meaningful  cathode  modelling  is  possible  with  this  formulation  as 
it  is  shown  in  reference  [13]. 

Applying  Eq.  (5.13),  the  energy  equation  (5.10)  is  transformed  into; 

f)  f)  f)  /  \ 

wND  +  aAWr.jr.)  =  +(j(r.) ,  (5,i5) 

where  coefficients  a,  b{Te),  c{Te)  and  Q{Te)  are  dehned  as 


Coefficients  5(Te)  and  c(Te)  depend  weakly  on  temperature  through  the  Hall  parameter  f3  = 
uje/ ^(.{Tfj).  Coefficient  Q{Te)  is  the  energy  source  term  and  accounts  for  (from  left  to  right  in 
its  dehnition)  cross-held  Joule  heating,  ionization  losses,  terms  related  to  density  variations 
along  the  streamlines,  and  energy  huxes  to  the  sheaths.  These  exchanges  with  the  sheaths  are 
obtained  from  the  wall  interaction  subcode. 
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Finally,  the  integration  of  equation  (5.9)  along  the  streamsurfaces  yields  a  (linear)  relation 
between  /e(A)  and  the  derivatives  of  Te(A)  and  0*(A), 


TedT^\ 
n*  dX  ) 


—  In 


(5.16) 


Equations  (5.14)-(5.16)  form  a  system  of  coupled  differential  equations  in  A  for  Te(A,t),  (j)*{X,t) 
and  /e(A,  t).  These  equations  depend  only  on  one  spatial  variable,  A,  thanks  to  the  assumptions 
of  constant  electron  temperature  and  Boltzmann  equilibrium  along  magnetic  held  lines.  The 
advantages  of  a  differential  formulation,  compared  to  the  integral  one,  is  that  it  gives  more 
hexibility  on  the  choice  of  the  most  convenient  numerical  integration  algorithm.  Also,  in  the 
new  code,  we  are  trying  to  make  the  magnetic  mesh  fully  independent  from  the  PIC-mesh.  In 
this  manner,  rehnement  can  be  achieved  in  the  electron  subcode  without  increasing  the  number 
of  cells  of  the  PIC-mesh.  Besides,  non-uniform  electron  meshes  can  be  used  to  better  model 
electron  physics.  Finally,  the  new  formulation  has  allowed  to  account  for  the  work  of  the  electric 
held  parallel  to  magnetic  held  lines,  neglected  in  HP  Hall-2. 

To  sum  up,  the  main  assumptions  and  features  of  the  macroscopic  electron  model  are  the 
following: 


1.  Te  is  constant  along  the  magnetic  streamlines,  so  that  Tg  =  Te(A)  with  A  dehning  the 
streamline. 


2.  A  generalized  Maxwell-Boltzmann  equilibrium  applies  for  the  electric  potential; 

e0(A,  a)  =  e(j)*{X)  +  Te(A)  In  (5.17) 

n* 

where  a  is  the  spatial  variable  along  a  streamline,  4>*{X)  is  the  ’thermalized’  potential,  n* 
is  a  reference  value,  and  ne(A,  a)  is  provided  by  the  particle  sub-code. 

3.  The  generalized  Ohm’s  law  relates  linearly  the  derivatives  of  Te(A)  and  4>*{X),  with  the 
electron  current  A  (A). 

4.  The  electron  current,  /e(A),  is  computed  from  a  current  conservation  ordinary  equation, 
which  admits  current  sources  from  the  walls; 


d 


(A  +  A) 


i=l,2 


j  •  n 
B-N  i' 


(5.18) 


where  A(A)  is  the  ion  current  function,  provided  by  the  particle  sub-code,  j  =  ji  +  je  is 
the  electric  current  density  at  the  inner(l)  or  outer(2)  walls,  n  is  the  unit  vector  pointing 
outwards  of  the  domain  and  N  is  dehned  in  Fig.  5.1. 

5.  The  temperature  prohle  Te(A)  is  obtained  from  a  parabolic  heat  equation  on  variables  A 
and  time  t.  This  equation  includes  A  (A)  and  heat  sources  ^q(A)  from  the  lateral  contour 
of  the  domain  (i.e.  lateral  thruster  walls  mainly). 
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5.1.4  Cross-field  electron  transport 

Classical  theory  predicts  an  electron  cross-field  transport  much  smaller  than  the  one  mea¬ 
sured  experimentally  in  HET.  Two  physical  mechanisms  have  been  proposed  to  explain  this 
high  electron  mobility;  anomalous  diffusion  due  to  correlated  azimuthal  ffuctuations  of  plasma 
density  and  electric  field[24],  and  near  wall  conductivity  due  to  secondary  electron  emission 
from  ceramic  walls[21,  25,  26].  These  two  mechanisms  can  be  incorporated  as  additional  elec¬ 
tron  cross-field  transport  through  effective  collision  frequencies.  Thus,  the  different  terms  in 
t'e  =  ^'en  +  +  t'ano  are  modelled  as  follows[3,  27]; 


8T, 


TTTTlp 


^ano  ^ano^e’) 


27rn  ■  je 


^ Him 


“27rnp 


^r(A) 


B 


dx 


(5.19) 


with  (Ten  =  3-10  and  aano  ~  10  ^  —  10  a  fitting  parameter  to  match  experiments  and 
simulation  results. 


5.1.5  Boundary  conditions  and  numerical  method 

In  order  to  integrate  electron  equations  proper  boundary  conditions  must  be  imposed  for  the 
electron  temperature,  the  thermalized  potential  and  the  discharge  current.  Electron  temper¬ 
ature  is  assumed  to  be  known  at  the  cathode  streamline  (5  eV)  while  a  zero  heat  conduction 
assumption  is  adopted  for  the  anode.  On  the  other  hand,  the  discharge  voltage  is  fixed  between 
the  anode  and  the  cathode,  and  consequently,  a  shooting  method  is  employed  to  verify  this 
condition  with  the  anode  potential  as  a  parameter.  Clearly,  the  anode  sheath  presented  in  a 
previous  chapter  must  be  accounted  for  in  these  computations[2]. 

The  numerical  method  we  employ  is  described  next.  Electron  equations  are  discretized  spa¬ 
tially  with  first-order  finite  differences  while  a  explicit  Euler  scheme  is  used  for  the  temporal 
discretization  (i.e.,  forward  in  time  and  centered  in  space,  ETCS).  Figure  (1.1)  shows  schemat¬ 
ically  the  magnetic  mesh  used  for  electron  equations  as  well  as  the  corresponding  magnetic 
field. 

PIC  sub-code  employs  a  time-step  At  such  that  most  particles  do  not  advance  more  than 
half  a  cell  per  time-step.  However,  electron  equations  must  be  integrated  with  a  smaller  time- 
step  St  due  to  numerical  stability  constraints  and  therefore  these  equations  must  be  solved 
several  times  per  particle  time-step.  In  fact,  the  ratio  At /St,  which  was  about  1.000  in  HPHall- 
2,  has  been  reduced  to  about  100  with  the  new  formulation.  This  change  has  more  than 
offset  some  increases  in  computational  cost  that  both  the  new  electron  formulation  and  the 
modifications  in  the  PIC-subcode  have  introduced.  So,  this  time  step  increase,  together  with 
other  improvements,  causes  a  20%  reduction  in  computational  time  and  permits  to  simulate 
1ms  of  thruster  operation  in  less  than  2  hours  with  a  2.5GHz  PC  computer. 
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5.2  Accuracy  improvement 


One  of  the  advantages  of  the  present  formulation  is  related  to  numerical  accuracy.  In  order 
to  evaluate  the  numerical  error  it  is  useful  to  exploit  the  electron  energy  balance  that  can  be 
expressed  as: 


P, 


elec.e 


+  Pc 


cathode 


=  Pr 


lateral, e 


+  Pc 


anode, e"! 


where  Peiec,e  is  the  overall  Joule  heating,  Pcathode  is  the  amount  of  energy  introduced  through 
cathode  in  the  form  of  an  electron  current,  Pioniz  represents  losses  due  to  ionization  and  excita¬ 
tion,  and  Piaterai,e  ^ud  Panode,e  represent  energy  losses  into  lateral  and  anode  sheaths  respectively. 
This  equation  is  derived  by  integrating  in  A  the  electron  energy  conservation  law  (5.10). 

Table  5.1  shows  the  electron  energy  balance  of  a  simulation  of  a  SPT-100  like  thruster.  The 
balance  is  almost  perfect  with  an  error  of  around  1%  of  the  total  electron  power. 


p 

^  elec,e 

P 

^  cathode 

P 

tomz 

P 

^  lateral, e 

P 

^  anode, e 

Error 

542.4  W 

12.8  w 

143.4  W 

326.5  W 

79.6  W 

5.7  W 

Table  5.1;  Electron  energy  balance  for  an  SPT-100  like  thruster 


It  is  possible  to  estimate  the  numerical  error  at  every  electron  time  step  based  on  electron 
energy  conservation  and  accounting  for  the  temporal  term  of  the  energy  equation.  The  evolution 
of  the  estimated  numerical  error  is  shown  in  hgure  5.3.  As  can  be  seen  it  almost  always  remains 
below  0.5  W.  Thus,  those  6W  of  error  in  the  stationary  energy  conservation  check  are  mainly 
due  to  the  temporal  term.  In  any  case,  it  is  around  1%  what  is  comparable  to  other  errors 
already  mentioned  in  this  report. 


Figure  5.3:  Time  evolution  of  the  estimated  numerical  error  associated  to  the  electron  energy 
equation 
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5.3  Near-anode  region 

The  simulated  spatial  domain  goes  from  the  external  plume  to  the  anode.  However  the  anode 
surface  A  (and  thus  the  anode  sheath  edge  B)  is  not  a  magnetic  streamline;  see  hgure  5.4.  Let 
S:  A  =  A5  be  the  magnetic  streamline  separating  those  streamlines  intersecting  the  anode  from 
those  without  intersection.  The  electron  model  is  not  well  designed  to  treat  the  near-anode 
region  between  S  and  the  anode  sheath  edge  B.  We  present  here  an  approximate  treatment  for 
this  region,  which  improves  clearly  previous  algorithms  used  by  HPHall-2. 

For  sake  of  simplicity,  here  we  consider  the  anode  (more  precisely,  the  anode  sheath  edge) 
as  a  z  =  const  surface.  The  new  algorithm  consists  on  the  following.  First  and  in  order  to  close 
the  spatial  domain  of  the  heat  equation,  we  consider  the  anode  as  a  magnetic  streamline.  The 
virtual  streamline  value  A*  assigned  to  the  anode  is  the  average  value  of  A(r)  along  the  anode. 
Second,  the  heat  and  Ohm’s  equations  are  solved,  thus  yielding  Te(A*)  and  0o(A*).  Third,  these 
values  are  not  assigned  to  the  anode  sheath  edge,  but  to  the  actual  streamline  A  =  A*,  which 
lies  between  B  and  S.  Four,  a  linear  interpolation  based  on  values  at  A5  and  A*  determines  Te{X) 
and  0o(-^)  in  the  whole  BS  region.  Finally,  Eq.(5.17)  yields  the  potential  0(A,(t),  wherefrom 
the  potential  distribution  at  the  sheath  edge,  (pni'f'),  is  obtained. 

The  new  algorithm  includes  improvements  in  the  model  for  the  anode  sheath  and  the  it¬ 
eration  procedure  to  determine  the  current- voltage  response,  IdiVd).  Let  (psir),  Tesir),  and 
nesi'f')  be  the  plasma  variables  at  the  sheath  edge,  obtained  from  the  integration  of  the  electron 
equations.  If  we  take  0  =  0  at  the  cathode  point,  the  potential  at  the  anode  would  coincide 
with  the  discharge  voltage  =  Vd- 


Figure  5.4:  Scheme  for  the  treatment  of  the  near-anode  region. 

The  model  of  the  anode  sheath,  presented  in  chapter  3  yields  a  relation  between  the  potential 
fall  (pAB  (at  a  given  location  r)  and  the  electron  current  jzeB  incident  onto  the  sheath, 

^(pAB  _  j  I 

TeB  \eneB\/TeB/mi 


(5.20) 
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The  algorithm  we  propose  considers  the  anode  potential,  0^,  as  a  known  parameter  and 
thus  the  anode  potential  fall  prohle  —  4>a  is  known.  This  algorithm  uses  the 

sheath  current- volt  age  curve  to  determine  the  anode  electron  density  current, 


jzeB{r)  =  eneB{r)‘ 


'TeBjr)  1  /  ecpABjr) 
mi  V  TeB{r) 


(5.21) 


Then,  in  order  to  satisfy  the  anode  sheath  conditions,  the  discharge  current  through  the 
external  circuit  should  be 


2TTdr[j^iB{r)  +  jzeB{r)], 


(5.22) 


where  the  ion  current  through  the  anode,  jziBi'f'),  is  provided  by  the  particle  code.  This 
current  R  should  coincide  with  the  current  injected  into  the  discharge  through  the  cathode 
(or  the  multiple  cathodes,  for  multi-stage  operation).  This  condition  is  used  to  compute  the 
anode  potential  <pA-  In  turn,  this  value  should  coincide  with  the  discharge  voltage  V^.  If  not, 
electron  equations  must  be  reintegrated  until  convergence  is  achieved;  2-3  iterations  normally 
are  enough. 


0.035  0.04  0.045  0.05  0.035  0.04  ,  ,  0.045  0.05 

r(m)  r{m) 

Figure  5.5;  Radial  prohles  of  plasma  magnitudes  at  the  anode  sheath.  CW,  SW,  and  BFC 
mean  corrected  weighting,  surface  weighting,  and  Bohm  forcing  condition.  Mj  is  the  ion  Mach 
number  at  the  sheath  edge.  Mj  and  Ue  are  computed  by  the  particle  code.  (pAB  is  obtained 
from  the  anode  sheath  model  of  Ref.  [2].  The  anode  potential  is  300V.  jeB  =  —je  ■  n 


At  present,  the  anode  sheath  model  is  valid  only  for  negative  or  electron- repelling  sheaths. 
It  could  happen  that  (pABi'f')  <  0  at  certain  zones  on  the  anode.  Provisionally,  for  these  zones, 
we  are  taking  from  Eq.(5.20)  the  value  of  jzeB  for  (pAB  =  0. 
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In  Fig.  5.5,  CW,  SW,  and  BCF  indicate  the  weighting  algorithms  used  by  the  particle 
code  to  weight  magnitudes  at  the  contour  nodes  [2].  The  Surface  Weighting  algorithm  seems 
the  best  and  more  consistent  option.  The  correct  value  of  ion  Mach  number  should  be  a  bit 
above  one.  The  decay  of  plasma  density  and  potential  fall  towards  the  lateral  wall  is  consistent 
with  the  gas  injector  being  at  the  center  of  the  anode.  Figure  5.6  shows  2D  maps  of  the 
near  anode  region  for  the  main  plasma  variables.  The  total  chamber  length  is  25  mm.  The 
separatrix  intersects  the  inner  and  outer  walls  at  6  mm  approximately.  For  this  case,  the  ion 
backstreaming  region  extends  only  ~  3  mm.  The  line  of  zero  ion  flux  coincides  approximately 
with  the  position  of  the  point  of  maximum  potential,  ~  319V.  The  large  radial  gradients  of  0 
for  BCF  in  Fig.  5.6  are  a  consequence  of  the  high  Mach  numbers  obtained  with  BCF,  which 
leads  to  excessively  low  plasma  densities  at  the  lateral  walls.  The  regular  variation  of  the  2D 
prohles  in  the  near-anode  region  can  be  presented  as  a  proof  of  the  suitability  the  algorithms 
presented  here.^ 
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Figure  5.6;  2D  maps  of  plasma  variables  in  the  near  anode  region.  The  left  contour  corresponds 
to  the  sheath  edge  B.  Sheaths  cannot  be  not  shown  in  this  (quasineutral)  scale. 


^One  interesting  observation  is  worth  here.  The  particle  sub-code  assumes  that  all  ions  arriving  to  a  sheath 
edge  are  collected  by  the  wall.  This  is  correct  for  negative  sheaths.  In  regions  with  positive  sheaths,  part  of  the 
ions  reaching  them  are  reflected  back  to  the  quasineutral  region.  Positive  sheaths  remain  to  be  implemented  in 
the  code.  It  should  not  be  diflicult  to  implement  in  the  PIC  subcode  the  reflection  of  ions  in  positive  sheaths. 
The  main  problem  lies  on  the  electron  model.  As  we  prove  in  Annex  A,  for  positive  sheaths  (and  even  for  small 
negative  sheaths)  inertial  effects  on  electrons  become  important  and  the  present  electron  formulation  is  unable 
to  include  them.  In  any  case,  positive  sheaths  are  not  common  in  normal  operation. 
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Chapter  6 

SIMULATIONS  AND  TESTS 


Recent  trends  in  Hall  thruster  design  have  shown  an  increasing  interest  on  dual-mode  efficient 
operation.  The  advantages  are  clear  for  this  type  of  operation;  on  the  one  hand,  a  high  thrust 
but  low  specihc  impulse  regime  allows  for  efficient  orbit  raising  where  the  minimization  of  the 
transfer  duration  is  the  main  constraint;  on  the  the  other  hand,  orbit  correction  tasks  such 
as  north-south  station  keeping  can  take  advantage  of  the  high  specihc  impulse  regime,  which 
minimizes  the  mass  consumption[28]. 

Several  programs  have  already  been  launched  in  the  last  few  years  to  design  and  qualify 
dual-mode  Hall  thrusters  and  the  main  problem  that  still  remains  is  the  operation  at  high 
specihc  impulse  (or  high  discharge  voltage).  It  is  usually  found  experimentally  that  a  decrease 
on  efficiency  appears  at  voltages  higher  than  600V-700V  and  it  is  claimed  that  one  of  the 
possible  reasons  is,  among  others,  the  excessive  generation  of  doubly-charged  ions,  leading  to 
ion  populations  of  very  diherent  average  velocities  at  the  thruster  exhaust. 

However,  recent  results  from  Hofer  and  Gallimore  [11]  have  shown  that  a  monotonically 
increasing  efficiency  is  possible  if  the  magnetic  held  topology  is  optimized  for  high  discharge 
voltages.  Thus,  disadvantages  of  multiple  ionization  seem  not  to  be  the  only  reason  for  the 
efficiency  decrease  at  high  specihc  impulse  operation  measured  in  other  experiments. 

A  promising  Hall  technology  for  dual-mode  operation  is  the  two-stage  design  where  an  in¬ 
termediate  electrode  is  used  to  enhance  the  ionization  and  acceleration  processes.  Although  the 
theoretical  advantages  of  this  type  of  thrusters  are  obvious,  experiments  tend  to  show  that  no 
clear  beneht  appears  when  two-stage  operation  is  used  [29,  30,  31].  This  might  be  due  to  a  poor 
understanding  of  the  physical  processes  involved  in  this  kind  of  thrusters  where  the  intermediate 
electrode  can  change  signihcantly  the  discharge  characteristics.  Ahedo  and  Parra[32]  published 
a  model  of  the  two-stage  discharge  where  they  demonstrate  that  performance  and  efficiency 
depend  much  on  the  position  and  current  of  the  intermediate  electrode  and,  technical  issues 
apart,  predict  relative  efficiency  enhancements  of  20%  for  optimal  values  of  the  intermediate 
electrode  parameters. 

This  chapter  is  aimed  to  analyze  the  results  the  new  code  is  capable  of  producing  precisely 
for  these  two  kinds  thrusters;  single-stage  high  specihc  impulse  and  two-stage  models.  Besides, 
some  parametric  variations  are  also  included  to  show  the  hexibility  of  the  simulation  code. 
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6 . 1  Performances 


6.1.1  Performance  characterization 

This  paragraph  is  intended  to  summarize  the  main  dehnitions  regarding  thruster  performance. 
When  we  analyze  the  efficiency  of  the  plasma  discharge  of  a  Hall  thruster,  the  natural  parameter 
is  the  so-called  anode  efficiency 


where  F  is  the  thrust,  Tjia  is  the  mass  flow  through  the  anode  and  is  the  electrical  power 
to  sustain  electrode  and  cathode  currents;  for  a  simple,  single-stage  thruster,  Pd  =  IdVd-  The 
anode  efficiency  can  be  split  up  into  several  partial  efficiencies  accounting  each  one  for  a  different 
process  inside  the  thruster.  Thus,  we  write 


(6.1) 


K  Tjy^  "fjcur  Vvol,i-\- 


(6.2) 


where 


Vu 


rrij 


m„ 


(6.3) 


is  the  mass  or  propellant  utilization  (subscripts  a  and  c  refer  to  the  anode  and  the  streamsurface 
containing  the  virtual  cathode,  respectively); 


^cur 


Idi,cCd 

Pd 


(6.4) 


is  the  current  utilization  (it  represents  the  additional  necessary  electron  current  to  maintain 
the  thruster  operation  and  can  also  be  considered  a  way  to  characterize  the  ionization  process); 


Pjet,i+ 

Vvol,i+  =  y  Ty 
J-di+,cVd 


(6.5) 


is  the  voltage  utilization  based  on  singly-charged  ion  acceleration  {Pjet,i+  is  the  singly-charged 
ion  jet  power  at  the  domain  exit); 


T;+/2ffi*+,c 

Vp,i+  =  —y - 


is  a  plume-divergence  efficiency,  based  on  singly-charged  ion  magnitudes  (Fj_|_  is  the  thrust  due 
exclusively  to  singly-charged  ions); 


(6.6) 


_  _ {^di+,c  A  Idi++,c  I^Y 

(^Idi+,c  T  Idi++,c)  (^Idi+,c  T  Idi++,c  /2)  ^  ^ 

is  the  charge  utilization  {Idi+,c  and  Idi++,c  are  singly-  and  doubly-charged  ion  currents),  which 
measures  the  influence  of  doubly-charged  ions  on  the  thrust  ;  and  K  accounts  for  minor  con¬ 
tributions  to  Tja  like  the  part  of  the  thrust  due  to  neutrals  and  electrons.  Sometimes  it  is  useful 
to  use  an  internal  efficiency  as  the  product  of  current  utilization  and  voltage  utilization. 
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6.1.2  Nominal  operation  results 

The  Hall  thruster  taken  as  reference  for  the  computations  is  of  the  SPT-100  type,  which  is 
characterized  by  the  nominal  properties  of  Table  1. 


Vd 

m 

Bfnax 

Pd 

7] 

F 

^sp 

300  V 

4.8  mg/s 

230  G 

1350  W 

0.48 

83  mN 

1600  s 

Table  6.1;  SPT-100  nominal  operating  and  performance  parameters. 


Pd 

7] 

F 

^sp 

1220  W 

0.39 

68  mN 

1434  s 

Table  6.2;  SPT-100  simulated  performance  parameters. 


The  simulation  results  obtained  with  the  conditions  sketched  in  Fig.  1.1  for  the  above  oper¬ 
ating  parameters  are  in  Table  2.  Many  other  averaged  performance  variables  can  be  found  in 
hgure  6.1  that  has  been  extracted  directly  from  the  performance  hie  generated  by  HPHall. 


Nominal  Voltage  [V] 

= 

300 

Thrust [N] 

=  0.0676132 

Voltage  [V] 

= 

300 

Current  [A] 

= 

4.06514 

Total  power  [W] 

=  1219.55 

Current  Oscillation  Level  [A] 

= 

1.01048 

Total  lateral  wall  loss[W] 

=  376.519 

Power  [W] 

= 

1219.55 

Thrust  power  [W] 

=  546.121 

Max  magnetic  field  (B-ave) [G] 

= 

249. 629 

Useful  power  [W] 

=  758 . 932 

Anode  magnetic  field (B-ave) [G] 

— 

33 .5621 

Thrust  efficiency 

=  0.389799 

Total  current  at  anode [A] 

= 

4.06514 

Internal  efficiency 

=  0.622304 

Anode  total  ion  current  [A] 

= 

-0 . 1687 

Specific  impulse [ s] 

=  1433. 9 

Anode  single  ion  current  [A] 

= 

-0.165812 

Anode  double  ion  current  [A] 

= 

-0.00288894 

PARTIAL  EFFICIENCIES 

Ion  current  to  lateral  wall [A] 

= 

1.57981 

Thrust  efficiency 

= 

0.389799 

Coefficient  Kl 

= 

1.08383 

Total  ion  beam  current[A] 

= 

3.08159 

Coefficient  K2 

= 

0. 992821 

Single  ion  beam  current[A] 

= 

2.77374 

Charge  efficiency 

0. 991879 

Double  ion  beam  current[A] 

= 

0 . 307846 

Plume  efficiency 

= 

0.711358 

Propellant  utilization 

= 

0.828559 

Nominal  Mass  flux[kg/s] 

= 

4 . 8e-006 

Current  utilization 

= 

0.758047 

Mass  flux  anode  line [kg/s] 

= 

4 . 80831e-006 

Voltage  utilization 

= 

0.817413 

Mass  flux  thruster  exit[kg/s] 

= 

4.81353e-006 

Beam  total  ion  mass  flux[kg/s] 
Beam  single  ion  mass  flux[kg/s] 
Beam  double  ion  mass  flux[kg/s] 
Maximum  ion  beam  mass  flux [A] 


3. 98397e-006 
3.77453e-D06 
2.09458e-007 
4.8e-006 


Figure  6.1;  Main  performance  variables  for  the  SPT-100  simulated  thruster 
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Clearly,  this  simulation  underestimates  thruster  efficiency.  There  are  two  kind  of  reasons. 
First,  no  attempt  was  made  to  optimize  the  ’empirical’  parameters  of  the  model,  like  the 
anomalous  diffusion  or  the  wall  accommodation  coefficients,  both  having  a  relevant  impact  on 
performances.  The  second  reason  lies  on  the  plasma-wall  interaction  model,  which  leads  to 
excessively  high  plasma  recombination  and  energy  losses  at  the  walls.  Depletion  of  the  the 
tail  of  primary  electrons[33,  34],  two  electron  temperature  models[26],  and  partial  trapping 
of  secondary  emission[35]  have  been  suggested  to  modify  the  common  plasma-wall  interaction 
model. 

Apart  from  this  overestimate  of  wall  losses,  the  code  is  capable  of  reproducing  well  the 
main  physical  aspects  of  the  discharge  and  can  be  used  to  analyze  (qualitatively,  at  least)  the 
influence  of  different  operating  parameters  as  discharge  voltage  or  mass  consumption. 

Figure  6.2  depicts  the  most  important  variables  for  the  nominal  operation  point.  These 
results  must  be  taken  as  a  reference  for  the  following  sections. 
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Figure  6.2:  Temperature  (Tg),  electron  density  (ue),  neutral  density  (u„)  (upper  hgures),  electric 
potential  (0(1/)),  ionization  rate  (hj)  and  ion  current  density  (jj^)  (lower  figures)  for  the  nominal 
operation  point. 
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6.1.3  Estimates  of  numerical  errors 

HPHall  is  a  non-stationary  code  with  a  strong  oscillatory  component  in  the  solntion.  Thus, 
all  ’time- averaged’  variables  have  an  associated  error.  This  analysis  is  intended  to  evaluate  the 
size  of  this  error  and  determine  the  best  trade-off  between  accuracy  and  computational  effort. 

To  this  end  several  runs  have  been  performed  with  the  same  operating  parameters  and 
similar  initial  conditions.  The  only  ’free’  parameter  in  this  study  is  the  number  of  iterations 
used  to  average  variables.  In  order  to  characterize  the  temporal  error  in  averaged  variables  we 
dehne  the  following  generic  variables:  the  time-average  and  the  typical  deviation. 
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Figure  6.3  shows  the  results  obtained  for  runs  with  10000,  20000  and  40000  iterations.  It 
can  be  seen  that  for  all  analyzed  variables,  the  relative  typical  deviation  is  always  smaller  that 
0.65%  what  implies  that  75%  of  simulations  have  an  uncertainty  smaller  than  1.3%.  This  value 
is  acceptable  and  consistent  with  the  numerical  errors  introduced  with  the  different  approxi¬ 
mations. 


Figure  6.3;  Relative  typical  deviation  of  different  variables  as  a  fnnction  of  the  nnmber  of  averaged 
PIC  iterations. 
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6.2  Parametric  studies 

In  this  paragraph  some  parametric  variations  are  presented  for  the  different  parameters  char¬ 
acterizing  thruster  operating  conditions.  To  this  end  we  start  from  the  SPT-100  case  analyzed 
in  the  previous  chapter  and  the  following  parameters  are  varied  continuously:  anomalous  diffu¬ 
sion  coefficient,  magnetic  held  strength,  mass  how,  thruster  channel  length.  Note  that  multiple 
ionization  is  not  considered  in  this  analysis  for  the  sake  of  simplicity. 


6.2.1  Anomalous  diffusion  coefficient 

In  order  to  ht  simulation  results  with  experimental  results  the  anomalous  dihusion  coefficient 
aano  ,  Eq-  (5.19),  is  normally  modihed  until  the  performance  is  optimized.  Figure  6.4  shows 
the  corresponding  results  for  this  parametric  variation.  Increasing  the  anomalous  dihusion 
coefficient  causes  the  discharge  current  to  increase  and  consequently  the  power  also  increases. 
Additionally,  the  higher  the  anomalous  dihusion  is,  the  lower  the  current  utilization  and  internal 
efficiency  are.  Thus,  the  overall  efficiency  decreases  dramatically  as  aano  increases. 


Partial  efficiencies 


Relative  specific  impulse  I  /I 

sp  SPj 


Figure  6.4;  Parametric  analysis.  Anomalous  dihusion  coefficient  variation  aano- 
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6.2.2  Magnetic  field  strength 

For  a  given  magnetic  circuit,  increasing  the  current  through  the  coils  causes  the  magnetic  held 
to  increase  proportionally.  Then  the  magnetic  held  gradient  also  increases  and  consequently 
the  thruster  operation  is  greatly  ahected. 

According  to  the  results,  increasing  the  magnetic  held  gradient  tends  to  reduce  the  accel¬ 
eration  region  and  to  move  it  downstream.  Thus,  it  is  closer  to  the  maximum  magnetic  held. 
Normally,  thrusters  are  designed  so  that  the  maximum  magnetic  held  is  at  the  thruster  chan¬ 
nel  exit.  That  is  why  an  improvement  in  the  overall  efficiency  is  observed  due  mainly  to  a 
better  current  utilization  and  lower  wall  losses.  However,  the  non-dimensional  specihc  impulse 

{Isp/Ispo  where  Rpo  =  ^  decreases  slightly  due  to  a  poorer  ionization.  As  expected, 

the  overall  power  decreases  because  a  lower  electron  current  is  required  as  the  magnetic  held 
strength  increases. 

Finally,  anode  losses  are  in  general  much  smaller  than  ceramic  wall  losses  since  electron 
temperature  near  the  anode  is  smaller  and  the  ion  back-streaming  how  is  much  smaller  the 
discharge  current.  Obviously,  it  is  also  important  to  take  into  account  that  the  anode  surface 
is  also  smaller  than  the  ceramic  wall  surface. 


Partial  efficiencies 


Efficiency  r] 


Relative  wall  losses  P  „  ,,/P , 
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Figure  6.5:  Parametric  analysis.  Magnetic  held  strength  variation. 
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6.2.3  Channel  length 

This  parametric  variation  consists  on  executing  the  program  for  several  thruster  geometries  with 
a  different  channel  length  keeping  the  rest  of  geometrical  parameters  constant.  The  channel 
length  is  one  of  the  most  important  variables  from  the  design  point  of  view.  As  it  can  be 
observed  in  hgure  6.6,  the  current  utilization  7]cur  decreases  and  the  propellant  utilization  7]u 
increases  due  to  the  relative  growth  of  the  acceleration  region  and  the  relative  reduction  of 
the  ionization  region  respectively.  Thus,  the  overall  efficiency  shows  a  peak  for  a  given  optimal 
channel  length  (around  30  mm  in  this  case).  As  in  previous  cases,  the  increase  of  the  propellant 
utilization  7]u  also  involves  higher  specihc  impulse  and  overall  power  (due  to  the  increase  of  the 
discharge  current  J^).  Finally,  as  expected,  wall  losses  increase  as  channel  length  increases. 
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Figure  6.6;  Parametric  analysis.  Channel  length  variation 
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6.2.4  Mass  flow  rate 

Figure  6.7  shows  the  variation  of  some  variables  characterizing  thruster  operation.  It  can  be 
observed  that  a  higher  mass  flow  rate  causes  an  increase  of  the  propellant  utilization  due  to  a 
more  efficient  ionization  enhanced  by  a  higher  neutral  density.  On  the  other  hand,  both  current 
utilization  and  internal  efficiency  decrease.  Because  of  this  behavior  of  the  partial  efficiencies 
the  overall  efficiency  shows  a  minimum  when  increasing  the  mass  flow  rate.  The  overall  power 
also  changes  due  to  the  increase  of  the  discharge  current,  which  is  proportional  to  the  mass 
flow  rate  and  the  propellant  utilization. 


m  (mg/s)  m  (mg/s)  m  (mg/s) 


m  (mg/s)  m  (mg/s)  m  (mg/s) 


Figure  6.7;  Parametric  analysis.  Mass  flow  rate  variation  m. 
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6.3  High  specific  impulse  thrusters 

This  section  presents  numerical  results  of  Hall  thrusters  operating  at  high-specihc  impulse.  The 
goal  of  these  simulations  is  two-fold.  On  the  one  hand,  we  aim  to  better  understand  the  plasma 
discharge  at  this  operation  conditions.  On  the  other  hand,  we  use  them  for  testing  the  new 
simulation  code.  Here  we  focus  on  results  for  high-specihc  impulse  and  single-stage  operation 
conditions. 

Figure  6.8  displays  the  variation  of  the  main  parameters  when  the  voltage  is  increased 
from  300V  to  750V,  and  the  magnetic  held  strength  is  increased  proportionally  to  the 
rest  of  parameters  remain  nominal.  As  expected,  the  specihc  impulse  and  the  thrust  increase 
proportionally  to  '  rather  well.  Deviations  are  due  to  variations  on  propellant  and  voltage 
utilizations.  The  total  discharge  current  remains  constant  through  the  entire  range  of  variation 
of  Vd  due  to  the  variation  of  the  magnetic  held  strength  commented  above.  As  a  consequence, 
the  total  power  varies  proportionally  to  the  discharge  voltage. 


150  300  450  600  750 


Figure  6.8;  Variation  of  main  operation  parameters  when  varying  the  discharge  potential  Vd-  Pwaii,i2 
and  Pwall,a  represent  total  energy  losses  to  inner  and  outer  walls  and  to  the  anode. 
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Efficiency  increases  up  to  600V  where  a  peak  is  found.  The  posterior  decrease  seems  to  be 
caused  both  by  higher  wall  losses  and  a  decrease  of  voltage  utilization  and  plume  divergence  effi¬ 
ciency.  The  same  parametric  investigation  has  been  performed  switching  off  multiple  ionization 
and  the  peak  of  efficiency  remains.  This  contradicts  the  common  argument  that  the  decrease 
of  efficiency  at  high  voltages  is  due  to  doubly-charged  ions,  but  seems  to  agree  with  Hofer  and 
Gallimore  who  obtained  a  monotonically  increasing  efficiency  by  optimizing  the  magnetic  held. 

Multiple  ionization  inhuence  is  quantihed  by  the  charge  utilization  efficiency,  77^,  and  the 
contribution  to  thrust  of  doubly-charged  ions  (Ti++/F).  At  750V  the  contribution  of  doubly- 
charged  ions  to  thrust  is  as  high  as  10%,  which  emphasizes  the  importance  of  including  and 
modelling  them  properly. 

Figures  6.9  and  6.10  display  prohles  at  the  channel  median  and  outer  wall,  respectively. 
Figure  6.11  shows  2D  plots  of  the  discharge  at  three  discharge  voltages.  In  all  cases,  time- 
averaged  magnitudes  (excluding  initial  transients)  are  plotted.  Relative  energy  losses  at  the 
walls  tend  to  decrease  when  increasing  Vd  except  for  the  limiting  case  of  750V.  This  reduction 
on  wall  losses  is  due  to  a  better  acceleration  process  which  prevents  the  ions  from  reaching 
the  walls  as  shows  the  prohle  of  the  ion  current  to  the  outer  wall,  jwi,  in  hgure  6.10.  The 
anomaly  encountered  at  750V  seems  to  be  caused  by  a  temperature  peak  inside  the  ionization 
chamber  as  it  is  shown  on  hgure  6.9.  This  peak  causes  higher  wall  losses  and  modihes  the 
electric  potential  fall  through  the  channel.  The  optimized  magnetic  held  topology  proposed  by 
Hofer-Gallimore  in  order  to  increase  efficiency  at  high  voltages  is  characterized  by  a  near  zero 
magnetic  held  strength  precisely  where  our  model  predicts  a  peak  on  temperature.  Thus,  it 
must  be  numerically  checked  whether  a  magnetic  held  topology  like  that  one  could  avoid  the 
peak  temperature  and  improve  efficiency  at  high  voltages.  Figure  6.9  also  plots  the  plasma 
density  and  equivalent  charge  number.  Due  to  electron  pressure  and  electric  acceleration,  the 
plasma  tends  to  be  conhned  near  the  anode  for  higher  voltages.  Thus,  the  ionization  process  is 
conhned  to  that  region  too,  as  shown  in  hgure  6.11. 

Figure  6.9  also  plots  the  plasma  density  and  equivalent  charge  number.  Due  to  electron 
pressure  and  electric  acceleration,  the  plasma  tends  to  be  conhned  near  the  anode  for  higher 
voltages.  Thus,  the  ionization  process  is  conhned  to  that  region  too,  as  shown  in  hgure  6.11. 
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Figure  6.9:  Electron  temperature  Tg,  electric  potential  (j),  electron  density  Ue  and  equivalent  charge 
number  Zi  profiles  at  different  discharge  voltages  300  V(solid),  600  V(dashed)  and  750  V(dashed- 
dotted)).  All  variables  are  evaluated  at  the  chamber  median  (r=0. 0425m).  The  vertical  dashed  line 
represents  the  thruster  exit. 
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Figure  6.10;  Electron  density,  Ue,  ion  current  to  wall,  jwij  energy  flux  to  wall,  qwi,  and  electron 
energy  flux  to  wall,  gvFe  for  Vd=  300  V,  600  V  and  750  V.  All  variables  are  evaluated  at  the  outer  wall 
(r=0.05m).  The  vertical  dashed  line  represents  the  thruster  exit. 
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Figure  6.11:  Electron  density,  ng  (upper  figures),  and  total  ionization  rate,  h*  (lower  figures)  at  = 
300  V,  600  V  and  750  V,  for  the  previous  parametric  variation. 
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6.4  Two-stage  Hall  thruster  simulations 

Figure  6.12  and  table  6.3  show  the  prohles  of  the  main  plasma  variables  and  performance 
parameters  when  an  intermediate  electrode  2mm  wide,  located  at  z=21mm  on  the  outer  wall, 
and  emits  different  electric  currents.  Rest  of  parameters  remain  nominal,  even  the  magnetic  held 
strength  and  prohle.  Both  propellant  and  voltage  utilization  tend  to  increase  when  augmenting 
the  electrode  current.  However,  since  the  magnetic  held  has  not  been  optimized  (as  Ahedo  and 
Parra  [32]  did)  the  current  utilization  and  thus  the  anode  efficiency  decrease. 
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Table  6.3:  Performance  parameters  for  diherent  currents  through  the  electrode  corresponding 
to  hgure  6.12. 


D 

h"i5 


no  electrode 

|\ 

0.26  A 

!  \ 

1.5E+18 

1 

0.79  A 

/  /  '  \ 

:  1.37  A 

w 

j:  /  '’W 

- 

7  lE+18 

i:/ 

5E+17 

. \ 

)  o.bl  0.02  0.03 

0.04 

% 

o.bl 

'y/\ 


300: 

250: 

^200: 

^150: 

100: 

50: 

0, 


8E+23|- 

6E+23 

A' 

'g  4E+23 
•  c" 

2E+23I- 


0  0.01  0.02  0.03  0.04 

z(m) 


”0  0.01  0.02  0.03  0.04 

z(m) 


z(m)  z(m) 

Figure  6.12;  Electron  temperature,  Te,  electron  density,  Ue,  electric  potential,  (/>,  and  total  ionization 
rate,  hi,  profiles  for  different  currents  through  the  electrode.  The  electrode  is  2mm  wide  and  is  located 
on  the  outer  wall  at  axial  position  z=21mm.  The  vertical  dashed  line  represents  the  thruster  exit. 
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Table  6.4:  Performance  parameters  for  diherent  electrode  locations  corresponding  to  hgure 
6.13. 

Table  6.4  displays  the  diherent  performance  variables  when  modifying  the  electrode  location 
and  maintaining  almost  constant  the  electrode  current.  It  is  clear  that  the  electrode  location 
ahects  the  thruster  response.  In  fact,  in  the  new  electrode  location  (z=16mm)  the  anode 
efficiency  is  higher  than  in  the  previous  case.  This  inhuence  agrees  well  with  the  conclusions 
of  Ahedo  and  Parra[32].  We  remind  again  that  the  current  utilization  is  very  low  due  to  the 
non-optimal  magnetic  held. 
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Figure  6.13:  Electron  density,  He,  ionization  rate,  rii,  and  electric  potential,  (j),  contours  for  reference 
case  (upper  figures)  and  a  case  with  an  electrode  2mm  wide  on  the  outer  wall  at  16mm  from  anode 
(lower  hgures). 
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6.5  An  external  cathode  model 

In  a  Hall  thruster  an  external  cathode  injects  the  electron  current,  Ic,  required  to  neutralize 
the  ion  beam  and  to  sustain  the  plasma  discharge  inside  the  thruster.  For  a  thruster  with 
ceramic  walls,  that  current  is  equal  to  the  current  through  the  anode,  and  it  is  known  as  the 
discharge  current.  Id. 

Let  C  be  the  external  boundary  (in  the  near  plume)  of  the  computational  domain,  which 
coincides  with  a  magnetic  streamline;  see  Fig.  6.14.  HPHall-2  considers  the  surface  C  as  the 
neutralization  surface.  First,  the  particle  sub-code  computes  the  ion  current  Re  crossing  that 
boundary.  Then,  the  electron  current  injected  inwards  of  the  domain  is  iRj  =  Iq  —  he-  This 
electron  current  is  a  boundary  condition  at  C  of  the  electron  equations  in  A.  A  second  boundary 
condition  is  the  ’cathode  reference  potential’  (0  =  0  in  our  plots,  generally)  set  at  a  point  on 
line  C.  A  third  boundary  condition,  required  by  the  heat  equation  is  the  temperature  at  C, 
Tec,  which  we  identify  as  the  temperature  of  the  injected  electrons.  HPHall-2  does  not  solve 
the  region  beyond  the  neutralization  surface. 

The  new  electron  model  presented  in  Ref.  [14]  admits  the  injection/collection  of  electrical 
current  through  the  lateral  walls  of  the  thruster,  illustrated  in  Eq.5.18.  This  enhancement  was 
envisaged  to  study  thrusters  with  two-stages  or  internal  segmented  electrodes,  but  it  is  clear 
that  it  can  also  be  used  for  a  more  flexible  and  realistic  treatment  of  the  electrode  injection 
process.  Hence,  we  decided  to  model  the  cathode  as  an  emissive  electrode  placed  in  the  external 
walls  of  the  thruster. 

Emissive  electrodes  require  a  different  sheath  model  for  this  type  of  electrodes.  This  is 
based  on  space-charge  saturated  sheaths  and  was  presented  in  Ref.  [36].  Figure  6.15  plots,  for 
an  emissive-electrode  sheath,  the  curves  relating  the  sheath  potential  fall,  (fwQi  fhe  net  electric 
current  density  jw  and  the  density  flux  of  electron  total  energy  through  the  sheath  edge  Q  (i.e. 
the  boundary  of  the  quasineutral  domain), 


Emissive 


Figure  6.14:  Scheme  for  the  new  cathode  model.  The  shadowed  region  is  the  neutralization 
layer.  The  potential  of  reference  is  set  at  the  wall  position  of  the  cathode. 
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jw/(eneQ(kVmJ^'")  jw/(eneQ(kT,/mJ^'=) 

Figure  6.15:  Emissive  electrode  response.  Relations  among  jw  =  j  -  n,  the  sheath  potential  fall 
(pWQ^  and  qIq  =  (f°Q  ■  n  for  ion  Mach  numbers  Mi  =  1,2,  and  3.  Electrons  are  emitted  cold  by 
the  emissive  cathode. 

The  neutralization  model  we  propose  is  the  following.  An  electrical  current  Ic  is  injected 
into  the  plasma  domain  through  an  annular  cathode  of  hnite  length,  placed  on  the  outer  external 
wall.  Fig. 6. 14.  The  cathode  surface  and  Ic  dehne  the  current  density  jjy.  The  sheath  model 
determines  the  energy  flux  qIq  and  the  sheath  potential  fall,  (pwQ-  The  fluxes  jw  and  q^°Q  are 
source  terms  for  the  electron  equations  in  A  (hence,  Ic  is  no  more  a  boundary  condition).  Setting 
the  ’cathode  reference  potential’  at  the  cathode,  the  potential  fall  (pwQ  provides  the  boundary 
condition  required  for  0  in  the  quasineutral  domain.  The  electric  current  I  at  the  boundary 
surface  C  is  now  set  to  zero  (i.e.  the  ion  beam  is  neutralized).  For  the  heat  equation  we  impose 
the  condition  dTf./d\\c  =  0,  which  we  intuit  more  correct  than  to  set  the  temperature  Te^. 

Figures  6.17  to  6.18  show  the  results  of  a  simulation  for  a  cathode  placed  between  r  =  60 
mm  and  r  =  64  mm  in  the  external  wall,  which  injects  about  ~  250mA/cm^,  resulting  in  a  total 
current  of  ~  4A.  The  magnetic  lines  intersecting  the  cathode  dehne  the  neutralization  layer. 
The  electrical  current  increases  from  zero  at  the  external  boundary  to  the  ~  4A  current  when 


z  (m)  z  (m)  z  (m) 

Figure  6.16;  Plasma  response  along  the  channel  median  for  the  case  of  Fig. 6. 17.  The  neutral¬ 
ization  layer  cuts  the  median  at  34. 6-41. 8mm.  (The  peak  on  I{z)  is  a  numerical  defect). 
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Figure  6.17;  2D  plasma  response  for  a  cathode  placed  in  the  external  wall.  Arrows  for  ji  and 
Un  do  not  scale  with  the  vectors  magnitudes. 


crossing  the  neutralization  layer  (Fig.  6.16).  Of  course,  this  variation  of  the  electric  current 
comes  mainly  from  the  injection  of  electrons,  R  changing  from  negative  to  positive  values  across 
the  neutralization  layer.  In  the  present  model,  we  impose  boundary  conditions  on  the  gradients 
of  Te(A),  but  the  value  of  Tg  is  imposed  nowhere  (in  the  quasineutral  plasma).  Hence,  we  hnd 
very  positive  that  the  simulation  yields  Tg  ~  8  —  lOeV  in  the  neutralization  layer,  in  good 
agreement  with  experimental  results.  That  energy  is  imparted  to  the  cathode  electrons  by  the 
potential  fall  at  the  cathode  sheath;  this  fall  is  about  10-15V  (Fig.6.18),  a  very  good  value  too. 
Along  the  channel  median,  plasma  potential  is  minimum  in  the  neutralization  region  (Fig.6.16), 
which  makes  sense,  since  injected  electrons  must  be  driven  out  of  that  region  (both  inwards  and 
outwards).  In  this  process  electrons  are  heated  by  Joule  effect.  The  strength  of  Joule  heating 
depends  on  the  plasma  perpendicular  resistivity,  which  is  proportional  to  the  magnetic  strength 
B  (or  to  if  Bohm  diffusion  is  negligible).  For  this  simulation,  the  magnetic  strength  is  200G, 
at  the  chamber  exit,  130G  at  the  cathode  position,  and  90G  at  the  external  boundary.  This 
would  explain  that  Joule  heating  is  similar  outwards  and  inwards  of  the  neutralization  layer. 
Glearly,  the  value  of  Tg  at  the  external  boundary  is  excessive.  We  expect  to  correct  this  defect 
by  placing  the  cathode  more  outwards. 

The  local  maximum  of  Ug  at  the  median  of  the  neutralization  layer  (Fig.  6.17)  would  be 
due  by  the  concentration  of  slow  ions  near  the  minimum  of  0.  Notice  in  Fig.  6.18  that  the 
plasma  density  at  the  thruster  external  walls  is  of  the  order  of  10^®m“^,  yielding  Debye  lengths 
of  0.33  mm,  thus  keeping  valid  there  the  hypothesis  of  quasineutrality. 
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Figure  6.18:  Plasma  response  along  the  inner  and  outer  walls  for  the  case  of  Fig.6.17.  U  refers 
to  the  location  along  the  wall.  The  sign  of  scalar  currents(j)  and  energy  fluxes  (g*"*)  to  the  walls 
is  determined  by  the  outwards  vector  n.  (pwQ  is  the  potential  fall  at  the  inner/outer  sheath. 
(pQMiz)  is  the  potential  fall  at  each  streamline  \{rM,  z)  between  the  channel  median  (r  =  vm), 
and  the  inner/outer  sheath  edge  Q. 


Appendix  A 

TRANSITION  FROM  NEGATIVE  TO 
POSITIVE  ANODE  SHEATH 


This  annex  has  been  submitted  for  publication  in  Physics  of  Plasmas  as  an  independent  pa¬ 
per  entitled  ’Presheath/sheath  model  of  the  Hall  thruster  near-anode  region  for  positive  and 
negative  falls’. 


A.l  Introduction 

There  exists  a  limited  understanding  of  the  plasma  behavior  in  the  near-anode  region  of  a  Hall 
thruster  discharge  and  its  influence  on  the  thruster  operation.  Experiments  and  models  have 
centered  the  research  on  the  ionization  and  acceleration  regions,  which  seem  more  relevant  for 
the  plasma  and  thruster  responses  and  are  more  accessible  to  direct  plasma  measurements. 
Nonetheless,  the  plasma  behavior  in  the  near-anode  region  is  important  in  determining  the 
anode  heating  and  can  affect  the  stability  of  the  whole  discharge.  Zhurin  et  al.  ,  in  a  review 
of  the  large  Russian  experience  on  Hall  thrusters,  assert  that,  under  normal  operation,  electron 
thermal  transport  is  more  than  sufficient  to  sustain  the  discharge  and  thus  a  negative  (i.e. 
electron-repelling)  sheath  needs  to  be  formed.  When  the  thermal  flux  is  insufficient  to  conduct 
the  discharge  current  a  positive  (i.e.  electron-attracting)  sheath  develops;  this  situation  arises 
for  low  propellant  flows  and  results  in  the  discharge  easily  becoming  extinguished 

Negative  sheaths  are  ion-attracting  and  result  in  a  noticeable  ion  backcurrent  in  the  rear 
part  of  the  chamber.  Ion  backcurrents  (existing  for  up  to  a  60%  of  the  channel  length)  were 
reported  experimentally  by  Bishaev  and  Kim  and  Kim  Direct  measurements  of  the 
near-anode  region  have  been  made  recently  by  Dorf  et  al  In  Ref.  [40],  positive  and 

negative  anode  falls  are  found,  depending  on  (i)  the  metallic  anode  being  clean  or  coated  by  a 
dielectric  him,  and  (ii)  the  values  of  the  control  parameters  (discharge  voltage  Vd,  mass  how, 
etcetera).  For  clean  anodes,  potential  falls  were  negative  preferentially  and  tended  to  decrease 
when  Vd  was  decreased.  In  Ref.  [41] ,  the  near-anode  region  for  clean  anodes  and  three  diherent 
magnetic  conhgurations  are  investigated.  Since  the  magnetic  lines  intersect  obliquely  the  anode, 
the  anode  potential  fall  is  not  uniform.  A  positive  sheath  was  measured  in  part  of  the  anode 
only  for  a  particular  magnetic  conhguration  (which  presented  an  intermediate  zero  magnetic- 
held  point).  Since  this  saddle  point  was  near  the  exit  of  the  chamber,  its  relation  with  the 
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behavior  of  the  anode  fall  is  unclear. 

The  fluid  and  hybrid(particle-in-cell/fluid)  models  of  Fife  Ahedo,  Martmez-Sanchez, 
and  co-workers[^]’[^^^’[^5,  and  Barral  et  al.  assume  the  existence  of  a  negative  sheath  and 
recover  the  ion  backcurrent  region  of  Bishaev-Kim  experiments.  In  agreement  with  the  obser¬ 
vations  of  Ref.  [40],  these  models  also  show  that,  when  Vd  decreases,  the  anode  fall  decreases 
and  eventually  vanishes.  However,  Ahedo  et  al.  pointed  out  that  electron  inertial  effects 
become  relevant  before  the  negative  sheath  vanishes,  thus  affecting  the  consistency  of  the  diffu¬ 
sive  approximation  generally  used  for  modelling  the  electron  fluid.  Ahedo  and  Rus^^^  proposed 
a  discharge  model  with  partial  electron  inertia  effects  and  demonstrated  that,  for  small  anode 
falls,  electron  inertia  limits  the  electron  azimuthal  energy  to  values  of  the  order  of  the  electron 
internal  energy,  thus  also  bounding  the  energy  deposition  at  the  anode.  Ahedo  and  Rus  used 
an  anode  sheath  model  proposed  by  Dorf  et  al.  ,  which  consists  of  two  regimes:  the  hrst  one 
corresponds  to  the  classical  large  negative  sheath  solution;  the  second  one  assumes  that  there 
is  no  anode  sheath.  The  transition  between  both  lies  at  the  point  where  the  electron  flux  from 
the  quasineutral  discharge  into  the  anode  sheath  equals  the  electron  thermal  flux  based  on  a 
Maxwellian  electron  velocity  distribution  function(VDF). 

However,  Dorf’s  no  sheath  model  presents  several  inconsistencies,  that  have  motivated  the 
present  work.  First,  the  assumption  that  the  electron  VDF  in  the  anode  sheath  approaches 
a  near-full  Maxwellian  distribution  is  valid  only  for  large  negative  sheaths.  As  the  negative 
potential  fall  tends  to  zero  the  electron  VDF  should  approach  a  half-Maxwellian  distribution, 
which  modihes  the  matching  parameters  of  the  sheath  and  presheath  solutions;  indeed,  it  will 
shown  that  a  complete  set  of  matching  conditions  between  the  partially  kinetic  sheath  and 
the  fluid  presheath  is  not  obvious.  Second,  there  is  no  clear  extension  of  Dorf’s  model  into 
a  positive  sheath  regime.  Third,  in  a  positive  sheath  a  kinetic  formulation  for  ions  must  be 
used,  whereas  the  electron  flux  perpendicular  to  the  sheath  is  expected  to  be  ’supersonic’  (in 
the  sense  given  by  the  appropriate  Bohm  condition),  which  means  that  electron  perpendicular 
inertia  terms  (not  included  in  Ref.  [8])  cannot  be  dismissed  in  the  anode  presheath. 

This  paper  presents  a  one-dimensional(lD),  two-region  model  of  the  near-anode  region, 
which  takes  into  account  full  inertia  terms  of  ions  and  electrons.  A  two-region  formalism,  based 
on  the  asymptotic  limit  ^ps,  with  A^  the  Debye  length  and  ipg  the  spatial  scale  of  the 

quasineutral  anode  presheath,  is  invoked.  The  goals  are  to  obtain  solutions  of  the  near-anode 
region  for  the  negative  and  positive  sheath  regimes,  and  to  demonstrate  the  existence  of  a 
continuous  parametric  transition  between  positive  and  negative  falls.  Section  H  discusses  a 
quasineutral,  fluid  model  for  the  anode  presheath  and  derives  a  generalized  Bohm  condition 
for  the  anode  transition.  Section  H  discusses  which  is  the  appropriate  kinetic  model  for  the 
repelled  species  in  a  weak  sheath.  Results  are  analyzed  in  Section  IV. 


A. 2  Anode  presheath 

Figure  A. 1(a)  sketches  the  two-region  model  of  the  near-anode  region.  The  anode  is  planar, 
metallic,  and  clean;  it  collects  all  ions  and  electrons  reaching  it;  secondary  particle  emission 
is  neglected.  The  anode  sheath  AB  is  seen  as  a  discontinuity  in  the  scale  ips  of  the  anode 
presheath.  Solutions  for  anode  sheath  and  presheath  must  be  matched  at  the  sheath  edge 
(point  B). 
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Figure  A.l;  (a)  Sketch  of  a  negative  anode  sheath  and  the  near-anode  region,  (b)  The 
cut-off  Maxwellian  of  (repelled)  electrons  at  the  negative  sheath  edge  B  {fza  =  0  for  Vza  < 
—  ^/2e(j)AB/'me)  and  at  the  anode  A  {f^a  =  0  for  Vza  <  0);  (j)AB  is  the  sheath  fall. 
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We  propose  a  presheath  model  that  omits  or  simplihes  those  aspects  and  phenomena  that 
are  not  essential  to  understand  the  near-anode  region.  Thus,  we  neglect;  (1)  plasma  production 
(by  either  ionization  or  wall  recombination);  (2)  variations  of  the  ion  and  electron  tempera¬ 
tures;  (3)  variations  on  the  electron  gyrofrequency  cjg  (>  0)  and  the  electron  (total)  collision 
frequency  z/gps-  Formally,  this  means  that  the  presheath  extension  £  is  smaller  than  the  spatial 
scales  characterizing  these  processes.  The  presheath  model  consists  of  equations  for  the  plasma 
density,  n  =  rif.  =  rii,  the  ion  axial  velocity  Uzi,  the  electron  (axial  and  azimuthal)  velocities, 
Uze  and  U0e,  and  the  electric  potential  0: 


^e^za  —  QzoiB  COUSt,  CX  ^7  65 
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Subscript  B  has  been  applied  to  magnitudes  that  are  assumed  constant.  All  these  constants 
(in  particular,  the  ion  and  electron  local  fluxes  and  the  local  Hall  parameter)  are  inputs  of  the 
presheath  model,  which  would  be  determined  from  the  matching  with  the  solution  for  the  rest 
of  the  thruster  chamber.  This  is  not  afforded  here;  thus,  the  presheath  equations  are  integrated 
from  a  location  Tar  from  the  anode’  (which  is  dehned  more  precisely  below)  to  the  sheath  edge. 
We  will  focus  the  analysis  on  situations  with  oJeB/i^eB  >  0(1)- 
Equation  (A.2)  yields  a  explicit  expression  for  (f){uzi), 


TiB  R^Uzi  =  const. 


(A.5) 


From  Eq.  (A.l)  and  quasineutrality  the  two  axial  velocities  are  related  by 


^22/^26  QziBldzeB' 


(AJ) 


In  the  Hall  thruster  discharge,  the  ion-to-electron  (i-e)  flux  ratio  can  be  expressed  as  Qzib! QzeB  = 
Ra/ Id  —  liAi  where  Id  is  the  discharge  (or  anode)  current  and  Ra  is  the  ion  back-current  at  the 
anode.  For  normal  thruster  operation  RA/Id  <  10%  typically,  and  gzw/gzeB  ~  Ra/ Id- 
Adding  Eqs.  (A. 2)  and  (A.3),  one  has 

[rrieul^  +  rriiuR  -  {T^b  +  Fj^)]  — =  me{uJeBU0e  -  VeBUze),  (A.7) 

^26  CtZ 

which  does  not  include  the  electric  held.  Equations  (A.4),  (A. 6),  and  (A.7)  determine  the  three 
relevant  velocity  components  of  the  plasma.  Notice  that  Eq.  (A.7)  does  not  include  the  electric 
force.  Thus,  as  long  as  oJeBl^eB  is  large,  the  effective  axial  accelerating  force  on  the  plasma  is 
the  magnetic  force  on  electrons. 

Since  Uze  cannot  become  zero,  Eqs.  (A.4)  and  (A.7)  state  that  the  only  singularity  (which 
we  identify  with  the  sheath  edge  B)  of  the  presheath  model  corresponds  to 

^^e^zeB  d~  ^^i^ziB  Hb  —  TeB  T  T^b- 


(A.8) 
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This  expression,  which  includes  both  UzeB  and  Uzib,  is  the  sonic  condition  for  this  plasma  with 
two  drifting  species.  If  both  sides  are  multiplied  by  UeB,  it  establishes  the  known  condition 
that  the  plasma  pressure  is  equal  to  the  plasma  axial  momentum  flux  at  a  sonic  point.  We  will 
refer  to  Eq.  (A. 8)  as  the  sonic  Bohm  condition  since  at  the  sheath  side  of  B,  Eq.  (A. 8)  must 
coincide  with  the  (marginal)  Bohm  condition  for  the  existence  of  a  valid  sheath  solution. 

It  is  important  to  observe  that,  were  d(j)/dz  an  external  electric  field,  Eqs.  (A. 2)  and  (A. 3) 
would  yield  that  the  ion  and  electron  fluxes  could  present  singular  points  at  Uzi  =  ^/TiB/nii 
and  Uze  =  \/TeB/me,  respectively.  However,  since  the  electric  held  is  self-adjusted  by  the 
quasineutral  plasma,  it  acts  as  an  additional  pressure  force  and  the  above  two  ’sonic’  points 
are  just  regular  points  of  the  plasma  how,  where  the  right-hand  side  of  Eqs.  (A. 2)  and  (A. 3) 
vanish.  In  particular,  Eq.  (A. 2)  states  that  the  electric  held  is  zero  at  Uzi  =  ^/TR^Jini. 

Using  Eq.  (A. 6)  and  the  sonic  condition  (A. 8),  the  ion  and  electron  axial  velocities  at  point 
B  satisfy 
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There  are  three  distinguished  regimes  for  the  near-anode  region; 
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Regime  N  corresponds  to  the  classical  low-speed,  dihusive-limit  for  the  electrons.  Regime 
P  represents  the  opposite  situation;  the  ion  velocity  (near  the  anode)  is  negligible,  plasma 
dynamics  are  dominated  by  the  electron  huid,  and  the  sonic  condition  determines  UzeB-  Regime 
I  represents  the  intermediate  case  with  ion  and  electron  kinetic  energies  of  the  same  order. 

Since  UzeUee  <  0,  the  right-hand-side  of  Eq.  (A. 7)  never  vanishes.  Therefore,  point  B  is 
necessarily  singular  with 
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These  derivatives  are  inhnite  in  the  quasineutral  scale  of  the  presheath.  Indeed,  they  are 
indicating  the  transition  to  the  much  smaller  spatial  scale  of  the  non-neutral  Debye  sheath, 
where  gradients  are  proportional  to  the  inverse  of  the  Debye  length.  Only  the  derivative  of 
U0e{z)  is  regular  at  point  B,  and  uoe  is  going  to  be  constant  inside  the  Debye  sheath. 

The  sign  of  d(l)/dz\B  determines  whether  the  anode  sheath  is  going  to  be  positive  or  nega¬ 
tive.  Equation  (A.  14)  states  that  there  is  a  single  no-sheath  case  separating  the  negative-  and 
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positive-  sheath  regimes,  which  is  characterized  by 

QziB  I 


QzeB 


T 


rriiT 


eB 


'^ziB  — 


iB 


T 


rrii 


^zeB  — 


eB 


nip 


For  the  no-sheath  case  the  electric  held  at  the  anode  is 
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that  is,  the  anode  attracts  electrons  mainly. 


(A.15) 


(A.16) 


A. 3  Anode  sheath 


According  to  the  above  presheath  model,  a  negative  sheath  is  expected  for 
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(A.17) 


and  a  positive  sheath  arises  for  the  opposite  inequality  signs.  Therefore,  for  both  sheath 
regimes,  the  attracted  species  enters  ’supersonically’  into  the  sheath,  whereas  the  repelled 
species  presents  a  ’subsonic’  drift.  Hereafter  we  present  the  formulation  for  a  negative  sheath. 
Because  of  the  symmetry  between  the  two  sheath  types,  the  solution  for  a  positive  sheath  comes 
out  from  just  exchanging  the  roles  of  ions  and  electrons.  Observe  that  the  symmetry  is  not 
perfect  because  of  the  electron  azimuthal  velocity,  U8e,  but  this  is  constant  within  the  sheath 
and  does  not  intervene  in  the  sheath  equations. 

The  thin  Debye  sheath  is  collisionless  and  unmagnetized.  Its  role  is  to  create  a  potential  fall 
that  adjusts  the  flux  of  the  repelled  species  that  is  collected  effectively  at  the  anode  to  the  flux 
that  comes  into  the  sheath  from  the  quasineutral  presheath.  A  kinetic  formulation  is  mandatory 
for  the  repelled  species.  Making  the  ansatz  that  the  sheath  potential  prohle  is  monotonic,  0  can 
be  used  as  the  spatial  variable  within  the  sheath.  Let  Vze  represent  the  velocity  of  individual 
electrons  perpendicular  to  the  anode,  and  fze{vze,<P)  characterize  the  ID  electron  VDF,  once 
integral  moments  on  the  parallel  velocity  components  have  already  been  carried  out. 

Since  electrons  satisfy  meV^^  +  2e0  =  const  in  the  sheath,  fzei^ze,  4>)  depends  only  on  the 
electron  VDF  at  the  sheath  edge  B,  fzeB{vze)-  This  distribution  must  satisfy  two  require¬ 
ments:  (1)  to  be  consistent  with  the  particle  behavior  in  sheath  and  anode,  and  (2)  to  match 
appropriately  at  the  sheath  edge  B  with  the  electron  fluid  formulation  in  the  presheath.  For 
instance,  the  usual  choice  of  fzeB{vze)  fully  Maxwellian  is  inaccurate  since  the  anode  does 

not  return  back  the  tail  of  electrons  reaching  it.  The  inaccuracy  is  certainly  negligible  when 
the  sheath  potential  fall  is  large  (i.e.  ciPab/Tsb  3>  1),  but  here  we  are  mainly  interested  in 
the  negative-to-positive  sheath  transition,  that  is  in  cases  with  ciPab/Teb  small.  The  drifted 
Maxwellian  (with  UzeB  as  drift  velocity)  proposed  in  Ref.  [6]  is  incorrect  too,  since  fzeB{vze) 
must  be  symmetric  on  Vze  for  electrons  that  are  reflected  back  within  the  sheath,  that  is  for 
I'^ze  I  <  A/2e0AB/uie- 

The  correct  choice,  within  the  ’family’  of  Maxwellian  distributions,  is  a  zero-drift,  cut-off 
distribution  [Fig.  A.  1(b)] 
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Here,  T*,  n*,  and  (p^B  are  the  three  parameters  characterizing  the  VDF,  and  H{vz)  is  the  Heav¬ 
iside  step  fnnction.  The  three  parameters  must  be  determined  from  three  matching  conditions 
at  point  B  between  the  plasma  formulations  in  sheath  and  presheath.  The  matching  conditions 
we  propose  are:  (1)  the  continuity  of  the  electron  density; 
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(2)  the  continuity  of  the  electron  flux; 


(A,20) 

and  (3)  the  fulhllment  of  the  marginal  Bohm  condition  at  the  sheath  edge  B; 

^[ne(0)  -  ni{<t))]B+  =  0+,  (A.21) 

with  ne(0)  and  ni(0)  the  density  prohles  in  the  sheath. 

The  general  Bohm  condition,  d{ne—ni)/d(f)B  >  0,  comes  out  from  imposing  that  the  solution 
of  the  Poisson  equation,  dP^cp/dz^  =  e(ne  — Uj),  departing  from  the  quasineutral  point  B  presents 
a  monotonic  potential  The  derivative  of 
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is  obtained  easily.  Then,  if  the  fluid  formulation  for  ions  is  kept  inside  the  sheath,  condition 
(A.21)  becomes 

neB  ^  exp{-e(pAB/T^)  _  n^B  23) 

Ti,  2T^^J^Te(pAB/T^,  ^i'^liB  ~  d^iB 

Thus:  Eq.  (A. 19)  relates  n*  and  UeBj  Eq.  (A. 20)  determines  (pAB  in  terms  of  QzeBi  and  Eqs.(A.23) 
and  (A. 8)  relate  and  T^b-  As  a  consequence,  the  plasma-wall  problem  can  be  solved  without 
having  to  determine  the  sheath  prohles. 

Let  US  check  that  the  present  model  recovers  the  classical  large  sheath  and  provides  a  good¬ 
matching  condition  at  the  positive-to-negative  sheath  transition.  For  ^  1,  Eq.  (A.  19) 

and  (A. 23)  yield  UeB  —  "n*  and  UziB  —  a/ (T*  -|-  Tib) /mi.  This  last  expression,  together  with 
Eq.  (A. 20)  [which  assures  that  UzeB  ^  A/T*/me]  and  Eq.  (A. 8),  yields  T^b  —  E*.  Therefore, 
for  e<pAB/TeB  ^  1,  fzeB  IS  close  to  a  full-Maxwellian  and 


UziB  — 


TeB  +  Ti 


iB 


nii 


UzeB 


(-tf)' 


(A.24) 


These  results  change  for  moderate  and  small  negative  sheaths.  In  particular,  for  e(pABlTi,  -C  1, 
fzeB  tends  to  be  a  half-Maxwellian,  UeB  —  'n*/2,  TeB  —  2T*/7r,  and  Uzw  and  UzeB  approach 
the  no-sheath  values  of  Eq.  (A. 15).  For  positive  sheaths,  Eqs.(A.18)-(A.24)  are  applicable  if 
subscripts  i  and  e  are  exchanged,  and  (pAB  is  substituted  by  —cpAB  =  I^abI- 
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Figure  A. 2;  Ion  and  electron  (axial)  velocities  versus  the  i-e  flux  ratio  at  the  sheath  edge.  The 
asterisk  represents  the  transition  from  positive  to  negative  sheath  for  Tis/TeB  =  0.2.  Vertical 
and  horizontal  dashed  lines  correspond  to  the  large  sheath  approximation.  The  dash-and-dot 
line  is  the  no-sheath  model  proposed  by  Dorf  et  al.  . 


The  continuity  of  the  sonic  Bohm  condition  across  the  sheath  edge  was  already  invoked  as 
a  matching  condition  in  the  presheath-sheath  model  of  Ref.  [10].  There,  different  macroscopic 
models  of  a  three-species  plasma  in  sheath  and  presheath  had  to  be  matched  at  the  sheath 
edge.  The  suitability  of  this  matching  condition  here  is  well  justified  too.  On  the  one  hand, 
it  is  based  on  the  continuity  of  a  physical  property,  the  ’sonic’  character  of  point  B  (i.e.  the 
equality  between  the  static  and  dynamic  pressures).  On  the  other  hand,  since  (f)AB  0  for 
TaB / TTia,  {oi  =  i,e),  it  recovers  exactly  the  full  ranges  of  the  positive  and  negative 
sheaths.  This  is  not  immediate  to  satisfy,  taking  into  account  the  different  plasma  formulations 
in  sheath  and  presheath;  the  models  of  Refs.  [42]  and  [6]  are  clear  examples.  Also,  in  a  previous 
version  of  this  model  we  imposed,  as  third  matching  condition,  that  the  electron  temperature 
obtained  from  the  VDF  [Eq.  (4)  of  Ref.  [2]]  should  coincide  with  the  electron  temperature  in  the 
presheath,  TeB-  R  turned  out  that  (j)AB  =  0  for  Uz^b! ^/TeB/^e  ~  0.9  so  that  no  satisfactory 
solution  could  be  found  for  the  range  0.9  <  Uzcb/ \/TeB/^e  <  1-  To  conclude,  it  can  be 
said  that  the  matching  of  the  integral  moments  of  fzeB  with  the  fluid  formulation  is  exact  up 
to  order  one,  and  only  partially  exact  in  order  two.  Therefore,  some  differences  at  point  B 
between  sheath  and  presheath  values  are  expected  in  high-order  plasma  magnitudes  that.  The 
good  news  are  that  for  the  main  macroscopic  magnitudes  the  Fortunately,  the  differences  are 
small  for  all  common  magnitudes.  For  instance,  the  flux  of  electron  energy  deposited  at  the 
anode  (a  relevant  magnitude  indeed)  is,  depending  on  the  sheath  regime. 
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One  has  that  Qzca/ QzeB  —  equal  to  nTeB  for  (pAB  — ^  0+  and  3Tes  for  Pab  — ^  0  , 

the  difference  being  a  mere  5%  at  the  positive-to-negative  sheath  transition. 
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Figure  A. 3;  (a)  Sheath  potential  drop  versus  the  perpendicular  macroscopic  velocity  of  the 
repelled  species  (called  a).  The  dashed  line  corresponds  to  the  standard  large  sheath  approxi¬ 
mation.  The  dash-and-dot  line  corresponds  to  the  model  proposed  in  Ref.  [2].  (b)-(c)  The  two 
constants  n*  and  T*,  dehning  the  cut-off  Maxwellian  at  the  sheath  edge,  fzaB- 
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Figure  A. 4;  Plasma  parameters  at  the  sheath  edge  versus  the  i-e  flux  ratio,  for  Tje/Tes  =  0.2 
and  oJesl^eB  =  10;  origin  of  0  placed  at  anode,  (pA  =  0;  Uq  =  gzeB{'^e/TB)~^^‘^ ■  The  dashed 
line  marks  the  no  sheath  case. 
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A. 4  Results 


Figure  A. 2  plots  the  relative  values  of  Uzi  and  Uze  at  point  B,  Eq.  (A. 9),  as  the  near-anode 
region  evolves  from  presenting  a  negative  to  a  positive  anode  fall.  The  dimensionless  curve 
of  Fig.  A. 2  is  just  the  universal  ’circle’  expressed  by  Eq.  (A. 8);  only  the  position  of  the  no¬ 
sheath  point  depends  on  the  ratio  Tj^/Tes  and  the  type  of  gas,  Eq.  (A. 15).  Dashed  lines 
represent  the  standard  ’large  sheath’  approximations,  Eq.  (A. 24),  for  positive  and  negative 
falls.  The  horizontal  dash-and-dot  line  is  the  no-sheath  regime  proposed  by  Dorf  et  al.  as 
continuation  of  the  (large)  negative  sheath  regime. 

Figure  A. 3(a)  shows  the  variation  of  the  sheath  potential  drop  with  the  drift  velocity  of 
the  repelled  species  from  the  ’large  sheath’  to  the  ’no  sheath’  cases.  It  can  be  observed  that 
the  ’large  sheath’  model  yields  already  an  error  of  a  20%  for  UzaB / a/Tqb /m^  ~  0.15.  The 
approximate  model  of  Ref.  [2]  is  plotted  too.  Figures  A. 3(b) (c)  depict  the  two  other  constants 
characterizing  the  cut-off  Maxwellian  at  point  B,  n*  and  T*  in  Eq.  (A. 18),  in  terms  of  UaB,  TaB, 
and  the  drift  velocity. 

Figures  A.4(a)-(d)  depict  the  evolution  of  plasma  magnitudes  at  the  sheath  edge  B  with  the 
i-e  flux  ratio  in  the  anode  presheath.  The  main  features  we  observe  are:  the  smooth  transition 
of  the  potential  fall  from  positive  to  negative  sheaths;  the  continuous  change  of  UziB  and  UzeB 
with  the  anode  fall;  the  decrease  of  UeB  as  we  move  towards  a  positive  sheath  (if  Qzeb  is  kept 
constant);  the  no  influence  of  the  local  Hall  parameter,  cuefi/zzes,  on  the  anode  fall;  and  the 
small  range  of  values  of  the  i-e  flux  ratio  leading  to  a  positive  sheath.  This  range  is  proportional 
to  the  square-roots  of  the  electron-to-ion  mass  ratio  and  the  ion-to-electron  local  temperature 
ratio,  Eq.  (A.  15).  In  a  real  Hall  thruster  discharge,  Tis/T^B  is  expected  small  but  Rb  cannot 
be  neglected  if  the  positive  sheath  regime  wants  to  be  recovered;  one  has  oc  Rb  for  a 

positive  sheath.  Also  the  ion  velocity  at  the  sheath  edge  can  never  be  taken  zero:  for  a  large 
positive  sheath,  the  small  ion  drift  plays  the  same  role  than  the  small  electron  drift  for  a  large 
negative  sheath. 

Figure  A. 5  depict  prohles  of  plasma  magnitudes  in  the  anode  presheath  for  different  values 
of  Szib/ QzeBi  covering  both  the  negative  and  positive  sheath  cases.  The  prohles  of  the  electric 
potential  present  a  maximum  within  the  presheath  where  the  ion  hux  satishes  Uzi  =  ^/TiB/rrii, 
Eq.  (A. 2).  As  \uziB  \  decreases  from  Rb / nii  (’large  sheath’  limit),  the  maximum  of  0(z)  moves 
toward  the  sheath  edge  B  and  disappears  eventually  in  the  positive  sheath  regime. 

Figure  A. 6  illustrates  the  behavior  of  the  electron  azimuthal  velocity.  It  plots  the  dependence 
of  U0eB  on  both  the  i-e  hux  ratio  and  the  local  Hall  parameter,  and  some  examples  of  spatial 
prohles  of  uoe  in  the  anode  presheath.  The  main  observation  is  that,  except  for  regime  N,  the 
azimuthal  velocity  dihers  substantially  from  the  dihusive  approximation 
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the  diherence  being  maximum  at  the  sheath  edge  and  increasing  as  oJ^B/^eB  increases.  Ahedo 
and  Rus  already  found  (within  a  model  of  the  full  Hall  thruster  discharge  with  partial  inertial 
ehects  on  electrons)  that  inertia  tends  to  limit  the  increase  of  the  azimuthal  energy  to  values 
of  the  order  of  the  thermal  energy.  The  behavior  of  uoe  in  regime  P  is  analyzed  below. 

The  limitation  of  the  azimuthal  electron  energy  are  good  news  for  the  heat  deposited  at  the 
anode  in  the  positive  sheath  regime.  Figure  A.  7  plots  the  energy  deposited  by  electrons  into 
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the  anode  for  negative  and  positive  sheaths  (the  energy  deposited  by  ions  is  much  smaller  in 
generally).  The  relative  heath  flux  Qzca/ OzesTeB  ranges  from  the  classical  value  of  2  for  regime 
N  to  ~  5  at  the  no-sheath  point.  Then,  in  regime  P,  since  the  azimuthal  energy  is  almost 
invariant,  the  changes  of  the  electron  energy  flux  are  due  to  \e(j)AB\,  which  is  proportional  to 
TiB- 


A. 4.1  On  the  character  and  extension  of  the  presheath 

The  presheath  prohles  of  Figs.  A. 5  and  A. 6  were  obtained  integrating  from  a  location  far  from 
the  anode,  where  ^  Tb  and  Eq.  (A. 26)  applies.  Whereas  the  solutions  for  neg¬ 

ative  and  positive  sheaths  are  equivalent  (except  for  n^e),  the  presheath  behavior  for  regimes 
N  and  P  presents  very  important  differences.  First,  there  is  the  small  range  of  parameters 
corresponding  to  regime  P  (mainly  for  heavy  gases).  Second  and  more  relevant,  the  presheath 
character,  which  is  determined  by  the  behavior  of  the  attracted  species  with  respect  to  the  ap¬ 
plied  magnetic  held,  is  totally  different:  regime  N  corresponds  to  an  ’unmagnetized’  presheath, 
while  regime  P  pertains  to  a  ’parallel-incidence,  magnetized’  presheath.  An  analytical  solution 
for  this  case  (and  UeB / ^eB  large)  was  derived  by  Ahedo  .  Appendix  A  recovers  that  solution 
with  some  improvements.  The  behavior  of  U0e{z)  in  regime  P  is  given  by  Eq.  (A. 32)  and  differs 
substantially  from  the  diffusive  behavior  for  regime  N,  Eq.  (A. 26).  From  Eqs.(A.32)  and  (A. 36), 
the  azimuthal  velocity  at  the  sheath  edge  is 
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This  asymptotic  solution,  plotted  in  Fig.  A. 6 (a),  shows  a  good  agreement  with  the  exact 
solution.  This  conhrms  its  suitability  to  interpret  the  plasma  response  in  regime  P.  For  instance, 
it  shows  that  the  anode  presheath  changes  from  being  collisional-diffusive  while  Uzej \/TB~f^e  < 
VeBl^eB  to  being  a  free-acceleration  zone  near  the  sheath,  when  Uzd =  0(1)-  Also, 
Eq.  (A.27)  conhrms  that  the  electron  azimuthal  energy  remains  of  the  order  of  the  thermal 
energy  except  for  \n{ujeB/veB)  1  (but  this  range  is  very  unlikely  in  Hall  thrusters,  at  least). 
Reference  [44]  yielded  only  the  dominant  (logarithmic)  term  of  the  right  hand  side  of  Eq.  (A.27); 
for  hi{ujeB/ ^bb)  not  very  large,  the  contributions  of  the  other  terms  is  not  negligible,  a  27%-40% 
in  Fig.  A. 6(a). 

The  extension  of  the  presheath,  also  differs  substantially  for  the  unmagnetized  presheath 
of  regime  N  and  the  parallel,  magnetized  presheath  of  regime  P.  An  immediate  observation, 
corroborated  by  Fig.  A. 5,  is  that  larger  forces  (i.e.  larger  gradients)  are  necessary  to  accelerate 
magnetically-trapped  electrons  in  regime  P  than  to  accelerate  unmagnetized  ions  in  regime  N. 
Figure  A. 8  plots  the  change  of  ips  and  the  presheath  potential  drop  A0ps  with  the  i-e  flux  ratio 
and  the  local  Hall  parameter;  these  values  are  based  on  the  criterion  that  the  anode  presheath 
extends  from  the  sonic  point  B  to  the  (arbitrary)  point  where  =  0.1Tb.  Observe 

that  A(j)ps  is  independent  of  the  Hall  parameter  and  has  different  signs  for  regimes  N  to  P, 
as  already  commented.  In  regime  P,  the  solution  of  the  Appendix  yields  that  the  presheath 
extension  scales  as 

(A.28) 
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Figure  A. 5;  Plasma  profiles  in  the  anode  presheath  for  different  values  of  QzibI 9zeB  =0.0001(1), 
0.00091(2),  0.01(3),  and  0.05(4).  Other  parameters:  RbIT^b  =  0.2  and  ijj^bIvcB  =  10.  Dashed 
lines  (3)  represent  the  no-sheath  case.  The  anode  sheath  is  seen  here  as  a  surface  discontinuity, 
black  points  indicating  values  at  the  sheath  edge  B. 
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Figure  A. 6;  (a)  Electron  azimuthal  velocity  at  the  sheath  edge  versus  the  i-e  flux  ratio  for 
^esI^eB  =10  and  30;  other  parameters  as  in  Fig.  A. 5.  The  dashed  line  marks  the  no  sheath 
case.  The  dashed-and-dot  lines  correspond  to  the  diffusive  approximation,  Eq.  (A. 26),  in  regime 
N,  and  to  Eq.  (A. 27)  in  regime  P.  (b)  Profile  of  the  electron  azimuthal  velocity  in  the  anode 
presheath  for  the  same  cases  than  in  Fig.  A. 5. 


Figure  A.  7;  Energy  deposited  by  electrons  at  the  anode  versus  the  i-e  flux  ratio,  for  for 
^esI^eB  =10  and  30;  other  parameters  as  in  Fig.  A. 5.  The  dashed  line  marks  the  no  sheath 


case. 
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Figure  A. 8;  Length  and  potential  drop  of  the  anode  presheath  versus  the  i-e  flux  ratio  for 
^esI^eB  =  10  and  30.  The  presheath  is  considered  to  extend  between  point  B  and  the  point 
where  —  0.1Tb. 


with  le  the  standard  electron  gyroradius.  Equation  (A. 28)  explains  the  small  dependence  with 
the  Hall  parameter  observed  in  Fig.  A. 8.  In  regime  N,  where  the  diffusive  limit  holds,  Eq.  (A. 7) 
becomes 


dUzi  ^  UjIb  +  vIb  meQzeB 

diZ  ^eB  ^^iQziB 


(A.29) 


and  the  presheath  extension  is 


^ps 


^^edzeB 


(A.30) 


which  scales  linearly  with  the  i-e  flux  ratio,  Fig.  A. 8(b). 


A. 5  Discussion 

We  have  presented  a  model  of  the  near-anode  region  of  a  Hall  thruster  that  shows  a  gentle, 
continuous,  parametric  transition  between  negative  and  positive  sheath  regimes.  A  central 
result  is  the  non-existence  of  a  no-sheath  regime  but  just  a  single  no-sheath  case.  Fundamental 
elements  of  the  model  are  (1)  to  include  full  inertia  terms  in  the  fluid  equations  of  ions  and 
electrons,  (2)  to  opt  for  a  cut-off  Maxwellian  distribution  of  the  repelled  species  in  the  sheath, 
and  (3)  the  use  of  the  (generalized)  Bohm  sonic  condition  as  the  last  matching  condition 
between  the  sheath  and  presheath  formulations.  Classical  solutions  for  large  negative  and 
positive  sheaths  (with  classical  Bohm  conditions)  are  recovered  asymptotically.  The  regime 
adopted  by  the  plasma  depends  on  the  relative  flux  of  ions  and  electrons  reaching  the  near¬ 
anode  region.  The  parametric  range  for  positive  sheaths  is  much  smaller  because  of  the  small 
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electron-to-ion  mass  ratio.  The  present  model  model  is  of  interest  to  any  device  or  problem 
where  large  variations  of  the  electron  and  ion  fluxes  are  likely. 

The  study  has  been  focused  in  situations  where  the  local  Hall  parameter  is  not  small,  and 
the  E  X  B  drift  (i.e.  azimuthal  drift  in  Hall  thrusters)  dominates  electron  dynamics.  Then, 
it  has  been  shown  that  the  diffusive  approximation  for  the  azimuthal-to-axial  electron  drift 
applies  only  in  the  large,  negative  sheath  regime.  For  the  positive  sheath  regime,  the  azimuthal 
electron  energy  remains  of  the  order  of  the  electron  thermal  energy,  which  limits  the  increase 
of  the  electron  energy  deposited  at  the  anode.  Also  in  the  positive  sheath  regime,  the  anode 
presheath  has  a  strongly-magnetized,  parallel-incidence  character,  which  reduces  greatly  the 
presheath  extension  (to  the  order  of  the  electron  gyroradius,  roughly).  The  case  of  a  small  Hall 
parameter  in  the  near-anode  region  has  not  been  discussed  in  detail  here,  since  it  is  simpler 
and  straightforward  from  the  present  analysis. 

The  scope  of  the  present  model  is  limited  to  the  near-anode  region.  Therefore,  it  can  teach 
nothing  about  how  the  main  parameters  governing  the  Hall  thruster  discharge  (i.e.  magnetic 
held  topology,  discharge  voltage,  mass  how  rate)  inhuence  the  plasma  response  in  the  near¬ 
anode  region.  In  particular,  the  coupling  of  this  region  with  the  rest  of  the  channel  would  yield 
how  the  operation  parameters  determine  the  i-e  hux  ratio  near  the  anode  and,  as  a  consequence, 
the  near-anode  regime.  From  existing  models  of  the  full  discharge,  strictly  valid  for  non-small 
negative  sheaths  only,  the  i-e  hux  ratio  is  known  to  decrease  when^^^’  the  discharge  voltage 
is  decreased,  or  the  magnetic  held  is  increased,  or  the  magnetic  held  shape  is  more  hat,  or  the 
turbulent  dihusion  decreases.  This  last  case  would  agree  with  the  assertion  of  Ref.  [41]  that 
the  turbulent  dihusion  is  lower  for  the  magnetic  held  conhguration  yielding  a  positive  anode 
fall. 

A  model  of  the  full  Hall  thruster  discharge  which  includes  all  electron  inertia  terms  (and, 
therefore,  the  negative-to-positive  transition)  faces  numerical  difficulties  that  have  not  been 
solved  yet,  even  in  the  one-dimensional  formulation.  Furthermore,  in  the  two-dimensional  case, 
the  region  near  the  anode  becomes  much  more  complex  since  huxes  are  not  uniform  and  the 
magnetic  held  impinges  obliquely  on  the  anode  and  connects  it  with  the  lateral  dielectric  walls. 


Appendix:  presheath  solution  for  large  positive  sheaths 


In  regime  P,  the  ion  drift  energy  can  be  dropped  from  Eq.  (A. 7),  yielding 


Tf, 


duz 


B  \  u-u-ze 

^ze  )  ;  ^eB^de  ^eB^zei 


nieUzeJ  dz 

which  is  solved  together  with  Eq.  (A. 4).  These  lead  to  the  hrst-integral 


Gb  (  ,  Tb 

~'^6e  I  ^ze  T 


^eB 


^^e^ze 


+ 


^Ib  + 


eB 


^eB 


{z  -  Zo) 


(A.31) 


(A.32) 


{zo  is  a  constant)  which  yields  uge  explicitly;  the  last  term  of  the  right-hand-side  is  going  to  dom¬ 
inate  for  oJeBl^eB  large.  Equation  (A.32)  is  to  be  compared  with  the  dihusive  approximation, 
Eq.  (A. 26),  for  regime  N. 

Substituting  Eq.  (A.32)  in  Eq.  (A.31),  one  has 


dz 


—  —  t'eB 


- 


(A.33) 
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where  the  dimensionless  variables  are 


u. 


'^ZP.  - 


^jTB/rrii 


t'eS  — 


^eB 


eB 


+ 


Z=  {z-  Zo)‘ 


eB 


'M^eB  +  ^cb) 
Tn 


(A.34) 


If  the  Hall  parameter  is  large,  i.e.  VeB  ^  1,  an  approximate  analytical  solution  of  Eq.  (A. 33) 
is  available  for  the  whole  range  0  <  Uze  <  1:  for  Uze  "C  1,  the  solution  is 


TT  2;"'  „ 

Uze[Z)  exp  yorfc 


and,  for  Uze  =  0(1),  the  solution  is 


rU  z^ 

^  -  Inu^e  +  y  -  -  ln(z>eBV^), 


(A.35) 


(A.36) 


where  the  constant  in  the  right  hand  side  has  been  determined  by  matching  Eqs.  (A.35)  and 
(A.36)  in  the  overlapping  region  z/es  Uze  ^  1,  where  both  solutions  are  valid. 

Setting  Uze  =  1,  Eq.  (A.36)  yields 


Zb  ^  -  ln(27r)  -  1, 

which  determines  zq  and  completes  the  solution. 


(A.37) 
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Appendix  B 

PLASMA- WALL  MODEL  FOR 
PARTIAL  THERMALIZATION 


This  annex  has  been  published  as  a  paper  entitled  ’Combined  effects  of  electron  partial  ther- 
malization  and  secondary  emission  in  Hall  thruster  discharges’  in  Physics  of  Plasmas  (vol.  14, 
083501,  2007). 

B.l  Introduction 

The  plasma  interaction  with  the  lateral  dielectric  walls  of  a  Hall  thruster  chamber  leads  to 
plasma  recombination,  deposition  of  ion  and  electron  energy,  and  deposition  of  electron  az¬ 
imuthal  momentum  (known  as  ’wall  collisionality’)  [37,  47,  21],  These  plasma  depositions 
affect  signihcantly  the  thruster  performances.  In  addition,  ion  impacts  on  the  walls  are  respon¬ 
sible  for  wall  erosion,  which,  at  the  end,  limits  the  thruster  lifetime.  An  accurate  model  of 
the  plasma  wall  interaction  is  one  of  the  main  open  problems  in  Hall  thruster  research.  The 
electron  velocity  distribution  function  (EVDF)  and  the  characteristics  of  the  secondary  electron 
emission  (SEE)  from  the  wall  are  two  key  aspects  of  this  problem. 

Hobbs  and  Wesson  [9],  assuming  a  Maxwellian  EVDF  and  cold  SEE,  showed  that,  as  the 
SEE  yield  increases,  the  potential  fall  at  the  Debye  sheaths  around  the  dielectric  walls  decreases 
and  the  electron  energy  deposition  increases.  For  a  SEE  yield  close  to  1  (~  98.3%  for  xenon), 
they  found  that  the  negative  (i.e.  electron-repelling)  sheath  reaches  the  charge-saturation 
limit  (CSL)  and  regime,  which  prevent  the  vanishing  of  the  negative  sheath  and  place  an  upper- 
bound  on  the  deposition  of  electron  energy  at  the  walls.  Ahedo  and  coworkers  applied  the 
Hobbs- Wesson  sheath  model  to  macroscopic  and  hybrid  two-dimensional  (2D)  models  of  the 
Hall  thruster  discharge,  and  found  out  that  it  yields  excessive  energy  losses  and,  consequently, 
unreasonably  strong  deterioration  of  thruster  performances  [21,  10,  3,  48]. 

Recent  works  suggest  that  the  plasma  is  not  collisional  enough  to  replenish  the  tail  of  high- 
energy  electrons  that  impact  the  wall,  so  that  the  distribution  function  of  bulk  electrons  is 
non-Maxwellian,  thus  presenting  a  smaller  ’effective’  temperature  in  the  direction  parallel  to 
the  magnetic  field  [49,  33,  26,  34,  50,  51].  Since  the  role  of  an  electron-repelling  sheath  next  to 
a  dielectric  wall  is  to  balance  the  electron  and  ion  fluxes,  a  partially-depleted  tail  of  the  EVDF 
reduces  certainly  the  sheath  potential  fall.  However,  the  effects  of  a  depleted  tail  on  the  ion 
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flux  itself  and  the  energy  losses,  the  two  main  magnitudes  of  plasma-wall  interaction,  need  to 
be  quantihed.  In  addition,  these  magnitudes  are  also  affected  by  the  amount  and  behavior  of 
the  SEE. 

Turning  now  to  the  characterization  of  the  SEE,  there  are  three  subjects  that  require  atten¬ 
tion;  (1)  a  correct  secondary  electron  emission  model;  (2)  the  trapping  and  thermalization  of 
the  SEE  beams  within  the  bulk  of  the  plasma;  and  (3)  the  magnetic  effects  on  the  SEE.  Both  the 
SEE  yield  and  the  energy  distribution  of  emitted  electrons  are  poorly  known  for  the  low  energy 
range  of  electrons  impacting  with  Hall  thruster  walls.  SEE  seems  to  be  constituted  by  slow  or 
’true-secondary’  electrons  (coming  from  internal  layers  of  the  wall)  plus  elastically  and  inelas- 
tically  reflected  electrons,  with  different  average  energies  and  emission  yields  [52,  53,  54,  55]. 
At  low  impacting  energies,  elastically  reflected  electrons  dominate,  with  an  emission  yield  of 
about  20-60%.  The  one-parameter  SEE  model  of  Hobbs- Wesson  considers  cold  true-secondary 
electrons  only.  In  a  model  of  the  full  discharge,  Barral  et  ah  [26]  took  into  account  SEE  with 
true-secondary  electrons  and  a  constant  yield  of  reflected  electrons  with  partial  energy  accom¬ 
modation;  no  parametric  investigation  was  carried  out.  Taccogna  et  ah  have  proposed  a  very 
detailed,  multi(>  10)-parameter  SEE  model  [55].  Here,  we  will  investigate  the  relevance  of  dif¬ 
ferent  parameters  characterizing  SEE,  in  order  to  propose  a  relative  simple  model  that  retains 
the  main  ones  only.  The  charge  saturation  limit  will  be  studied  too. 

Fife[16]  and  Ahedo[10]  assumed  total  thermalization  of  the  SEE  beams  within  the  main 
plasma.  This  justihed  the  use  of  a  single,  Maxwellian  electron  population  in  the  quasineutral 
plasma,  even  for  high  SEE.  Later,  Ahedo  and  Parra  [35]  pointed  out  that  weak  collisionality 
could  allow  to  partial  recollection  of  the  SEE  beams  by  the  walls,  and  demonstrated  that,  when 
the  recollected  fraction  is  large,  the  sheaths  do  not  reach  charge-saturation  and  electron  energy 
losses  to  the  walls  are  much  lower.  Indeed  the  losses  are  similar  to  the  zero  SEE  limit,  i.e. 
about  100  times  lower  than  for  total  thermalization  and  a  charge-saturated  sheath.  Sydorenko 
et  al.[51,  56]  recover  this  behavior  of  SEE  at  low  thermalization  with  a  particle-based  model. 

Magnetic  effects  on  true-secondary  electrons  have  been  considered  partially  by  Sydorenko 
et  al.  too.  They  included  the  effect  of  the  E  x  B  (azimuthal)  drift  imparted  to  beam  electrons 
once  they  leave  the  (thin)  emission  sheath,  and  pointed  out  that,  in  the  weakly-collisional  case, 
this  increment  of  energy  facilitates  the  partial  recollection  of  SEE  and  increases  the  re-emission 
yield. 

The  goals  of  this  paper  are  (1)  to  derive  an  quasi-analytical  model  of  plasma- wall  inter¬ 
action  that  takes  into  account  the  combined  effects  of  partial  depletion  of  the  main  electron 
population  and  partial  recollection  of  the  SEE  beams,  and  (2)  to  investigate  which  are  the 
main  SEE  characteristics  that  must  be  retained  in  a  SEE  analytical  model.  Emphasis  will  be 
put  on  (a)  comparing  the  low  and  high  thermalization  limits,  (b)  evaluating  the  different  roles 
of  secondary  and  primary  electrons  in  the  response,  (c)  the  sensitivity  of  CSL  conditions  to 
different  parameters,  and  (d)  deriving  simple  expressions  for  particle  and  energy  losses,  which 
can  be  implemented  in  full  models  of  the  plasma  discharge,  such  as  HPHall[3,  16].  Progresses 
on  this  work  were  presented  in  two  conference  papers  [57,  58]. 

As  in  similar  works,  the  analysis  of  the  collisional  processes  determining  the  state  of  the 
EVDF  remains  out  of  this  paper  scope.  The  approximate  expressions  used  for  the  EVDF 
within  the  sheaths  will  be  based  on  hrst-principles  models.  Magnetic  effects  on  the  SEE  will 
not  be  included  in  the  present  model.  Strictly,  this  limits  the  model  validity  to  (i)  magnetic 
lines  perpendicular  to  the  wall  and  (ii)  a  weak  enough  axial  (i.e  parallel  to  the  wall)  electric 
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field,  Ez-  Although  there  are  uncertainties  on  the  value  of  the  azimuthal  velocity,  uq  EJB, 
(which  is  measured  indirectly  only)  experiments  and  simulations  would  suggest  that,  in  most 
Hall  thrusters,  meUg  -C  3Te,  except  in  a  very  localized  zone  at  the  center  of  the  acceleration 
region,  where  Ez  is  maximum. 

B.2  Formulation  of  the  electron  model 

A  stationary  plasma,  confined  between  two  planar,  ceramic  walls  is  considered  (Fig.B.l).  A 
one-dimensional  model  is  proposed  (which  corresponds  to  the  radial  direction  in  a  Hall  thruster 
chamber);  when  necessary,  axial  plasma  contributions  are  included  as  source  terms.  The  zero 
Debye  length  limit  is  invoked,  leading  to  a  two-scale  structure,  consisting  of  the  bulk  region  of 
quasineutral  plasma  and  two  collisionless,  space-charge  sheaths  adjacent  to  the  walls.  Let  point 
M  be  the  channel  median,  points  W  and  W’  the  two  walls  and  points  Q  and  Q’  the  edges  of  the 
two  sheaths.  The  sheaths  are  semi-infinite  regions  in  their  natural  scale  (the  Debye  length)  and 
discontinuity  surfaces  in  the  quasineutral  scale  (so  that  tq  =  rw  =  h/2  and  tq/  =  rw'  =  —h/2). 
We  consider  only  sheaths  with  monotonic  potentials,  which  means  that  the  model  will  not  go 
beyond  the  CSL.  Then,  the  electric  potential  0,  instead  of  r,  can  be  used  as  the  independent 
variable  inside  the  sheath  (thus  avoiding  to  define  a  different  spatial  variable  for  the  sheath). 
Let  (pwQ  =  0Q  ~  (pw  be  the  potential  drop  in  the  sheath,  which  is  part  of  the  solution  if 
the  walls  are  dielectric.  In  order  to  neglect  all  magnetic  effects  in  the  model,  we  assume  that 
the  magnetic  field  B  is  radial  and  near-uniform,  and  the  azimuthal  drift  velocity,  ue,  is  much 
smaller  than  the  electron  thermal  velocity,  Ce  =  ^/T^/mP. 

The  velocity  of  each  electron  is  divided  into  components  parallel  and  perpendicular  to  1^, 
V  =  Vrlr  +  and  the  electron  distribution  function  has  the  functional  form  f{r,Vr,v±). 

The  symmetry  of  the  problem  with  respect  to  the  median  M  implies  that 

/(rM  +  Ar,  Vr,  v±)  =  /(rju  -  Ar,  -Vr,  v±).  (B.l) 


B.2.1  The  electron  distribution  function 


The  EVDF  /(r,  v)  is  obtained  from  the  one-dimensional  (radial)  Boltzmann  equation 


V  ^  + 

dr  rrie  dr  dvr 


(B.2) 


where  /  accounts  for  collisional  processes  (e-e  thermalization,  e-n  elastic  collisions,  ionization), 
the  effect  of  plasma  instabilities,  and  transverse  (i.e.  axial)  plasma  diffusion  (required  anyway 
to  balance  the  net  production  or  loss  of  plasma). 

In  the  zero  Debye  length  limit,  /  can  be  dropped  within  the  two  thin  sheaths  and  the  EVDF 
depends  only  of  the  constants  of  motion.  These  consist  of  the  the  parallel  and  perpendicular 
components  of  the  electron  total  energy; 

nj_  =  const,  men^/2  —  60  =  const.  (B.3) 


Therefore,  the  EVDF  in  the  sheath  has  the  functional  form 


(B.4) 


f{(f>,Vr,V±)  =  fQ{VrQ{(l>,Vr),V±)  =  /w (t'rW (0 Wr) ,  W) , 
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(a) 

W  Q'  M  Q  W 


Figure  B.l:  (a)  Sketch  of  the  radial  model  with  the  bulk  region  and  the  two  Debye  sheaths, 
(b)  Electron  populations  in  the  sheath  and  secondary  electron  emission. 

where  0  is  used  as  independent  variable, 

VrQ(4>,Vr)  =  sign(n^)y^n2  +  2e(0Q  -  0)/me,  (B.5) 

and  Vrw{.(l>Wr)  satishes  a  similar  expression.  Then,  we  just  need  to  determine  /q(i’)  or  fwiv) 
to  know  the  EVDF  in  the  sheath. 

A  consistent  determination  of  the  EVDF  inside  the  sheath  requires  solving  the  Boltzmann 
equation  in  the  bulk  quasineutral  region.  Even  using  the  simple  Bhatnagar-Gross-Krook  for¬ 
mulation  [59]  for  /,  the  problem  faces  two  big  difficulties:  hrst,  the  quasineutrality  condition 
implies  that  the  potential  prohle  0(r)  is  part  of  the  solution,  and  second,  /  depends  on  velocity 
moments  of  /(r,  v).  Based  on  phenomenological  considerations  and  the  simple  analysis  carried 
out  in  Ref.  [57],  we  assume  that  the  distribution  function  at  the  sheath  edge  Q  can  be  expressed 
as 

{ffQ{Vr,V±),  Vr<-VWQ, 

fliv),  \vt\<vwq,  (B.6) 

ftQ{Vr,V±),  Vr>VwQ, 
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where  vi^q  =  ^J2e(l)wQ/Tni,  ftq  corresponds  to  electrons  impinging  into  the  wall  W,  fjq  cor¬ 
responds  to  electrons  emitted  or  reflected  at  wall  W,  and  /i  corresponds  to  the  population 
of  primary  (or  bulk  or  thermalized)  electrons.  The  population  //  of  electrons  from  the  wall 
depends  on  the  SEE  properties  of  wall  W  and  will  be  discussed  below.  A  Maxwellian  function 
is  used  to  define  the  primary  population, 


=  ni 


312 

exp 


mev‘^\ 


(B.7) 


where  ni  and  Ti  are  constants  (notice  that  the  actual  distribution  function  of  primary  electrons 
is  determined  below).  The  population  at  Q  of  electrons  going  to  the  wall  satisfies 


/ie(«r,t>±)  =  (1  -<T)ffQ'(v,;V_L)  +  afi(v), 


(B.8) 


where  the  first  term  in  the  right-hand  side  is  the  non-thermalized  fraction  at  Q  of  electrons 
coming  from  wall  W’,  the  second  term  corresponds  to  primary  electrons,  and 


a{v)  ^  1  -  ex.p[-h/\ther{v)[ 


(B.9) 


is  the  function  representing  the  thermalization  fraction,  that  depends  on  the  thermalization 
mean-free-path,  Xther{v),  and  the  channel  width  h.  In  the  absence  of  a  reliable  model  for  the 
electron  thermalization,  we  take  ^{v)  =  const  as  in  Ref.  [51]. 

Secondary  electron  emission  at  wall  W  depends  on  the  distribution  of  electrons  impacting 
the  wall, 

ftw{Vr,V±)  =  ae~^^^fi{v)H{vr)  +  (1  -  (^)ffW'{Vr,V±)-,  (B.IO) 

here,  Eqs.(B.4)  and  (B.8)  have  been  used,  (fwQ  =  ^4>wq/Ti,  and  H{vr)  is  the  Heaviside  step 
function.  For  the  purposes  of  the  present  work,  only  the  main  aspects  of  the  experimental 
data  on  SEE  are  modelled.  Thus,  the  SEE  yield  is  assumed  to  have  contributions  of  elastically 
reflected  electrons  and  ’true-secondary’  or  beam  electrons. 


5s{E)  =  +  5sr{E),  (B.ll) 

with  E  =  mn^/2  the  energy  of  the  impacting  electron.  Subscripts  r  and  b  stand  for  ’reflected’ 
and  ’beam’,  with  the  last  name  justified  on  their  emission  energy  (~  1  —  3  eV)  being  generally 
much  less  than  the  sheath  potential  fall  and  the  temperature  of  primary  electrons  (~  5  —  50 
eV).  Based  on  experimental  data  [52,  53,  60,  61]  and  previous  models  [26,  55,  62],  the  partial 
SEE  yields  are  assumed  to  follow  the  simple  linear  laws 


Ssb  =  E/Eb,  Ssr  =  Soexp{-E/Er)]  (B.12) 

in  particular,  the  exponentially-decaying  law  for  6sr{E)  is  suggested  by  results  from  Refs.  [52, 
53,  55].  A  constant  value  for  Sgr  was  used  by  Barral  et  ah  [26].  A  possible  small  threshold  energy 
(of  few  eV)  in  Sgb  has  been  neglected  because  of  its  small  impact  on  the  response.  Typical  values 
for  Boron  Nitride  ceramics,  used  in  Hall  thrusters,  would  be  ~  Ef,  ~  50eV  and  (5o  ~  0.4  — 0.6, 
Eb  ~  40eV.  Notice  that  the  cross-over  energy  for  the  total  SEE  yield,  Ei,  comes  out  from 


(B.13) 
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The  distribution  of  electrons  from  wall  W  is  modelled  as 

ffw{Vr,V±)  =  Ssr{E)ftw{-Vr,V±)  +  f2{v)H{-Vr), 
where  the  hrst  term  in  the  right-hand  side  corresponds  to  reflected  electrons  and 

„  ,  ,  f2TmiI  (  rUe  /  mev‘^\ 


(B^14) 


exp(^-  — 

models  the  distribution  of  emitted  beam  electrons,  with 

g2=  [  d^VVrSsb{E)[l  -  5sr{E)]ftw{Vr,V±) 


(B.15) 


(B.16) 


representing  their  particle  flux.  Obviously,  only  electrons  collected  by  the  wall  can  produce 
true-secondary  emission.  The  emission  temperature  T2(~  1  —  3eV)  is  the  fourth  parameter  of 
the  SEE  model,  and,  as  commented  before,  T2/T1  is  small  in  most  or  all  the  thruster  chamber. 

The  populations  of  electrons  to  and  from  the  wall,  Eqs.(B.lO)  and  (B.14),  can  now  be 
expressed  as  linear  combinations  of  distributions  fi  and  /2, 


fiw(vr,v±)  =  [<xe  +  (1  - 


(B.17) 

(B.18) 


(B.19) 


1  -  (1  -  a)dsr 

is  a  gain  factor  representing  the  cumulative  effect  of  the  subsequent  collisions  of  reflected  elec¬ 
trons  with  the  walls,  Fig.B.2.  The  two  contributions  to  ftw  in  Eq.(B.17)  are  the  replenished 
tail  of  primary  electrons,  and  the  total  number  of  beam  electrons.  Since  one  has  aK  <  1  for 
a  <  1  and  ^Sgr,  partial  thermalization  always  implies  a  partial  depletion  of  the  incident  tail  of 
primary  electrons.  On  the  contrary,  (1  —  a)K  can  be  larger  than  one  (in  particular,  for  a  -C  1 
and  Ssr  >  0). 

We  can  now  express  the  electron  distribution  function  at  any  location  inside  the  sheath 
(dehned  by  its  local  potential  0)  as  the  sum  of  the  contributions  of  the  distributions  of  primary 
and  beam  electrons. 


=  /p(0,'n)  +  /6(0,'n). 

(B.20) 

/p(0WrW±)  =  X  <1  1, 

I'^rl 

(B.21) 

[  CTK, 

Vr  ^ 

{ 

0 

1 

V 

/s(0Wr,w)  =  X  0, 

\Vr  1 

(B.22) 

[  (l-a)K, 

Vr  ^ 

where  Vro  =  ^/2e{(j)  —  (j)w)/'fTie 


B.2.  FORMULATION  OF  THE  ELECTRON  MODEL 
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Figure  B.2:  (a)  Evolution  of  the  distribution  function  in  the  sheaths  for  a  =  0.6,  T2/T1  =  0.1, 
6sr  =  0.25,  Eb/Ti  =  1.5,  and  ^wq  =  1.91.  Dimensionless  variables  are  r  =  r/h,  v  =  Vr \/me/Ti, 
and  f{r,Vr)  =  ^jTi/nie  dv±27iv±f{r,Vr,v±).  Regions  <  0  and  >  0  correspond  to 
distributions  //  and  ft,  respectively,  of  electrons  from  wall  and  to  wall,  (b)  Comparison  of  the 
EVDF  at  the  two  sheath  ends,  Q  and  W.  (c)  Illustration  of  the  evolution  of  the  distribution 
function  in  the  bulk  region  (for  0  =  const) 


B.2. 2  Fluxes  of  particles  and  energy 

The  fluxes  of  particles  and  energy  at  the  wall  satisfy 

9e  =  J  d^VVrfw{Vr,V±), 

(B.23) 

QeW  =  J  d^VVr]^mv‘^fw{Vr,VA_). 

(B.24) 

Notice  that  particle  fluxes  are  constant  across  the  sheath  but  energy  fluxes  are  not;  at  the 
sheath  edge:  qeQ  =  QeW  +  (pWQQe- 

Because  of  the  presence  of  6sr  in  function  k,  the  above  integrals  cannot  be  expressed  in 
a  simple  non-integral  form  except  for  6sr{E)  =  const  (i.e.  E  —  r  =  00)  or  a  =  1.  We  solve 
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next  the  problem  for  6sr  =  const  and  postpone  the  more  general  case  to  Sec.  IV.B.  Using  Eqs. 
(B.17)  and  (B.18),  one  has 


9e  9p  9hi  QeW  2Ti5'p  2,T2gbi 


(B.25) 


with 

=  (1  -  5sr)(JKe~^^^ni\/Ti/2Time,  9b  =  0-^92,  (B.26) 

the  net  fluxes  or  primary  and  beam  electrons,  respectively.  Notice  that  g2,  Eq.(B.16),  is  only 
the  flux  of  ’new’  beam  electrons;  g2  and  gb  coincides  only  for  total  thermalization,  <7  =  1. 

For  Ssb{E)  linear,  this  flux  of  emitted  beam  electrons,  satishes 

92  =  —  Ssr)9tw/Eb,  (B.27) 


where 


9p  2Ti  2T2 

qtw  =  z - ^ — —  +  K{l-a)—g2, 

1  —  Osr  Eb  Eb 


(B.28) 


is  the  energy  flux  incident  into  the  wall.  Solving  Eq.(B.27)  for  g2,  the  net  beam  (or  ’true¬ 
secondary’)  yield  (dehned  as  the  ratio  between  beam  and  primary  net  fluxes)  is 


9p 


72p 


2Ti 


Eb  -  k{1  -  a)(l  -  6sr)‘2T2 


(B.29) 


where  72^  =  g2/ 9p  is  the  average  emission  yield  of  true-secondary  electrons. 

Appendix  A  gives  expressions  for  other  magnitudes  of  interest,  such  as  the  partial  densities 
of  primary  and  beam  electrons,  Up  and  Ub-,  respectively,  and  the  temperature,  T^pg,  of  primary 
electrons,  in  the  direction  parallel  to  the  magnetic  held  (it  is  immediate  that  the  perpendicular 
temperature  of  primary  electrons  is  Ti). 


B.3  Closure  of  the  sheath  model 


The  parameters  involved  in  the  EVDF  model  can  be  divided  in  several  groups.  The  hrst  one 
consists  of  the  type  of  gas,  which  dehnes  the  ratio  ^Jmi/nie-,  and  the  reference  values  ni,  Ti. 
These  are  used  to  dehne  non-dimensional  variables: 


—  ^ 
ni  ^  m^jTi/mi 


f 


T 


(B.30) 


and  so  on.  The  second  group  consists  of  the  thermalization  factor  a  and  the  four  parameters 
of  the  SEE  model;  T2,  Ef,,  Er,  and  (5o.  Finally,  there  is  the  sheath  potential  fall  (pwQ- 

If  the  potential  fall  is  known  (for  a  conducting  wall,  for  instance)  the  model  of  Sec.  II 
would  be  complete  in  order  to  determine  the  EVDF  and  the  electron- wall  interaction.  Figures 
B.2(a)-(b)  plot  an  example  of  the  evolution  of  the  EVDF  inside  a  sheath,  with  0  acting  as 
spatial  variable.  Figure  B.2(b)  compares  the  EVDF  at  the  two  sheath  ends;  beam  electrons 
and  the  partially  depleted  tail  of  primary  electrons  are  clearly  observable  in  fg.  Figure  B.2(c) 
has  illustrative  purposes  only:  it  shows  the  prohles  of  the  EVDF  in  the  bulk  of  the  plasma 
for  constant  electric  potential  there.  Observe  the  two  counterstreaming  beams  of  secondary 
electrons  and  their  partial  thermalization. 


B.3.  CLOSURE  OF  THE  SHEATH  MODEL 
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For  dielectric  walls,  the  potential  fall  (pwQ  is  determined  from  the  condition  of  zero  electric 
current  at  the  wall, 

9e  =  9i=  jiw/e,  (B.31) 

with  jiw  the  ion  density  current  into  sheath  and  wall.  This  condition  couples  the  ion  and 
electron  problems.  Since  jiw  is  going  to  depend  weakly  on  a  and  the  SEE  model,  it  is  convenient 
to  rewrite  the  expressions  of  Eq.  (B.25)  in  the  form 


9p 


9i 

1  -7fep’ 


2Ti  —  2T27fep 

QeW  =  — ^ - 9i- 

1  -  Ihp 


(B.32) 


Substituting  gp  in  the  hrst  equation  with  Eq.(B.26)  and  solving  for  the  sheath  potential  one 
has 

exp0^Q  =  —(1  -  -  7.p).  (B.33) 

V  27rme  9i 

The  ion  current  at  the  sheath  edge  depends  on  the  plasma  behavior  in  the  bulk  region.  As 
an  example  and  in  order  to  close  the  model,  we  determine  next  jiw  for  the  simplest  case  of  a 
stationary  regime  and  singly-charged,  quasi-cold  ions. 


B.3.1  A  simple  ion  model 

In  a  stationary  regime,  the  ion  current  into  the  sheath  is  determined  self-consistently  from  the 
Bohm  condition  at  the  sheath  edge.  The  Poisson  equation  for  the  sheath  potential  is 

^  —  K(0)  -  ni(0)] ,  (B.34) 

dC  eo 

where  the  ion  and  electron  densities  depend  only  on  0.  It  is  well  known  that  the  development 
of  a  space-charge,  monotonic  solution  from  point  Q"*"  (on  the  sheath  side)  requires  to  satisfy 
the  generic  Bohm  condition 

^  [ne(0)  -  ni(0)]Q+  >  0.  (B.35) 

Furthermore,  if  the  response  is  stationary  in  the  quasineutral  region,  only  the  marginal  (or 
sonic)  form  of  the  Bohm  condition  applies. 

The  density  of  a  quasi-cold,  singly-charged  ion  population  satishes 

ni{(l))  =  9i  [{9i/neQf  +  2e(0Q  -  (t))/m^]  ,  (B.36) 

where  quasineutrality  at  Q  has  been  applied.  Then,  the  marginal  case  of  Eq.(B.35)  leads  to 


9i  ^eQ'^riQi 


(B.37) 


which  is  further  developed  in  Eq.(B.49)  of  the  Appendix.  This  expression  states  that,  for  a  cold 
ion  population,  gi  depends  only  on  the  electron  density  around  Q.  For  the  simple  case  of  zero 
SEE  and  a  non-depleted  Maxwellian  distribution,  the  Bohm  condition  states  that  ions  enter 
the  sheath  with  the  sound  velocity  ^jT^Jmi. 
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B.3.2  The  charge  saturation  limit 

The  above  plasma  model  assumes  that  the  electric  held  is  monotonic  in  each  half  channel. 
This  condition  is  no  longer  valid  when  the  electric  held  in  the  sheath  becomes  zero  at  the  wall 
boundary,  i.e.  d(j)/dr\\Y  =  0,  which  is  known  at  the  charge  saturation  limit  [9].  The  integration 
the  Poisson  equation  across  the  sheath  shows  that  the  condition  d(j)/dr\w  <  0  is  equivalent  to 

r4>Q 

/  e(nj  —  ne)d(l)  >  0, 

J  (f>W 

that  is  the  net  electric  charge  inside  a  monotonic  sheath  is  non-negative.  This  equation,  detailed 
in  the  Appendix,  will  be  used  to  determine  the  CSL.  The  plasma  response  beyond  the  CSL  [54] 
is  not  considered  in  this  paper. 


(B.38) 


B.4  Results 

The  main  magnitudes  of  plasma-wall  interaction  affecting  Hall  thruster  performances  are:  (i) 
the  plasma  flux  to  the  wall,  Qi  (=  Qe),  Eq.(B.32),  which  is  recombined  and  requires  to  be 
reionized;  (ii)  the  deposition  of  electron  energy,  qew,  which  is  the  main  source  of  plasma-cooling 
and  wall-heating;  and  (iii)  the  ion  energy  flux  into  the  wall, 

QiW  =  {miU^iQ/2  +  ecpwQ)  Qi,  (B.39) 

which,  apart  from  producing  additional  wall  heating,  is  the  responsible  for  wall  sputtering.  The 
present  radial  model  gives  information  only  on  the  flux  of  radial  energy  qiw,r  =  +  c^iwq)  9 

B.4.1  Influence  of  the  thermalization  level 

Figure  B.3  depicts,  for  the  simplest  SEE  model  {Sgr  =  0  and  T2/T1  1),  the  combined 

influence  of  the  electron  thermalization  level,  a,  and  SEE  cross-over  energy,  Ei{=  Ei,),  on  the 
main  plasma-wall  variables.  Both  parameters  affect  strongly  the  plasma  response.  The  plasma 
flux  Qi,  governed  by  Bohm  condition,  remains  close  to  one  always  (except  for  a  — 0).  The 
sheath  potential  fall,  (pwqio'),  presents  a  maximum  for  an  intermediate  thermalization  level. 
This  behavior  of  (pwQ  is  easy  to  understand  from  Eq.(B.33);  The  hrst  factor  in  the  right-hand- 
side  of  Eq.(B.33)  is  large  (~  200  for  xenon),  the  second  one  is  the  effect  of  partial  thermalization 
on  the  primary  electron  population,  and  the  last  one  is  the  contribution  of  secondary,  beam 
electrons.  Therefore,  for  any  of  the  last  two  factors  to  affect  the  value  of  (pwQ  signihcantly,  they 
must  be  small.  This  occurs  for  either  n  1,  i.e.  the  low-thermalization  limit,  or  1  —  'fbp  ^  1, 
which  corresponds  to  the  vicinity  of  the  charge-saturation  limit.  The  density  of  beam  electrons 
at  the  sheath  edge  (and  in  the  bulk  region)  is  negligible.  A  partial  exception  is  at  the  CSL, 
where  UbQ  is  the  responsible  for  the  increment  of  Qi. 

The  most  dramatic  dependence  with  a  corresponds  to  the  deposition  of  electron  energy,  qew, 
which  increases  from  q^\Y/2Tigi  ~  1  at  a  -C  1  to  qew/‘^Tigi  ~  100  at  the  CSL.  This  increment 
is  due  to  the  increase  of  the  flux  of  primary  electrons  gp,  Eq.(B.32).  The  evolution  of  qiw,r{o') 
follows  that  of  (j)wQ-i  with  a  maximum  at  an  intermediate  thermalization  level  (which  can  lead 
to  maximum  wall  sputtering  there). 


B.4.  RESULTS 
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In  our  first  approach  to  this  model  [57],  we  assumed  a  constant-frequency  model  (i.e.Xther/'Ure  = 
const)  instead  of  a  constant-mean-free-path  one.  It  can  be  checked  now  that  both  models  lead  to 
the  same  trends  in  the  solution,  with  the  present  model  working  with  much  simpler  expressions. 


Figure  B.3;  Influence  on  the  plasma-wall  response  of  the  electron  thermalization  level  and  SEE 
cross-over  energy  for  Eb/Ti  =  1.5,  2,  2.5,  and  oo,  T2/T1  =  0.01,  and  6sr  =  0  (thus  Ei  =  Ef,). 
Asterisks  correspond  to  the  CSL.  This  and  following  hgures  are  for  xenon  ~  490). 


B.4.2  Influence  of  the  emission  model 

Figure  B.4  shows  the  effect  of  a  hnite  temperature  of  true-secondary  electrons,  T2.  The  increase 
of  T2/Eb  yields  a  larger  Ssb  and  thus  a  larger  'jbp-  This  leads  to  larger  gp  and  qew,  and  lower 
(j)WQ-  The  exception  is  a  ~  1,  when  the  only  effect  of  increasing  T2  is  to  reduce  Qew-  The 
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parameter  T2/ E^  is  small  always  and  we  expect  solutions  for  T2/£'fe  ~  0  to  yield  errors  of  order 
0{T2/Ei,).  In  Sec.  VI  we  will  see  that  this  is  no  longer  true  for  the  CSL,  where  errors  are  of 
order  0{^%fE^). 


Figure  B.4;  SEE  model.  Effect  of  the  emission  energy  of  ’true-secondary’  electrons:  T2/E11  = 
0.01,  0.1,  and  0.3;  Eh/Ti  =  1.8  and  6sr  =  0.  Asterisks  correspond  to  the  CSL. 


The  presence  of  elastically  reflected  electrons  modihes  the  sheath  parameters  more  deeply, 
since  these  electrons  are  hot  and  6sr  is  not  small  always.  Figures  B.5(a)-(d)  show  the  effect  of 
a  SEE  yield  with  a  0-60%  fraction  of  elastically  reflected  electrons.  Notice  that  the  plots  differ 
on  the  SEE  energy  that  is  kept  constant:  Et  in  Figs.  5(a)-5(b)  and  Ei  =  [1  —  6sr)Eb  in  Figs. 
5(c)-5(d).  The  different  plasma  behavior  is  due  to  d'jbp/dSsr  being  positive  for  Eb  =  const, 
and  negative  for  Ei  =  const.  This  illustrates  that  in  order  to  understand  correctly  the  plasma 
response  we  must  interpret  hrst  the  role  and  relevance  of  the  several  plasma  conditions  and 
parameters.  The  analysis  shows  that  the  zero-electric-current  and  Bohni  conditions  are  those 
affecting  the  most  the  plasma  response.  Since  Qi  varies  weakly,  the  key  parameter  (for  T2/Eb 
small)  is  the  relative  beam-to-priniary  net  electron  flux,  'jbp  =  Qp/ Qbi  which,  through  Eq.(B.32), 
self-determines  the  primary  flux  gp  and  therefore  the  energy  losses  Qew-  The  sheath  potential 
fall,  Eq.(B.33,  is  then  self-adjusted  in  order  the  wall  collects  the  appropriate  flux  of  primary 
electrons. 

Up  to  here  we  have  assumed  6sr{E)  =  const  in  order  to  obtain  simple  relationships  among 
the  parameters  that  facilitate  their  interpretation.  There  is  no  difficulty  in  obtaining  results 
with  Er  hnite  in  Eq.(B.12):  expressions  for  Qp  and  'ybp  just  keep  their  integral  form.  Figure  B.6 
plots  some  solutions  with  Er  hnite.  As  expected  there  are  no  qualitative  changes  with  respect 
to  the  case  6sr  =  const.  This  leads  us  to  propose  an  approximate  emission  model  that  keeps 
the  dependence  of  6sr  on  the  energy  of  impacting  electrons  but  proht  simultaneously  from  the 
simple  expressions  of  Secs.  II  and  III.  The  idea  is  to  substitute  the  ’exact’  yield  Ssr{E)  by  an 
’average’  one,  6sr,  which  depends  directly  on  the  main  electron  temperature  Ti.  A  satisfactory 
enough  agreement  is  reached  using  the  approximate  SEE  yield 

(B.40) 

in  the  analytical  expressions  of  Secs.  II  and  III.  The  dashed  lines  of  Fig. B.6  correspond  to  this 
approximate  model. 


B.  5.  THE  LOW- THERMALIZATION  LIMITS 
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Figure  B.5;  SEE  model.  Effect  of  the  presence  of  elastically-reflected  electrons:  6sr  =  0,  0.3 
and  0.6;  T2/T1  =  0.01;  and  E^/Ti  =  1.8  in  (a)-(b),  Ei/Ti  =  1.8  in  (c)-(d).  Asterisks  correspond 
to  the  CSL. 


B.5  The  low-thermalization  limits 


For  total  thermalization,  the  present  model  recovers  the  results  of  Ref.  [10].  The  case  of  low 
thermalization  is  of  special  interest  here  because  (i)  it  presents  important  differences  with  total 
thermalization,  (ii)  it  admits  some  asymptotic  expressions;  and  (iii)  it  could  be  the  appropriate 
one  for  Hall-thruster  discharges  (a  value  of  a  ~  1%  is  suggested  in  Ref.  [56],  which,  for  Ti  =  30eV 
and  a  chamber  width  of  15mm,  corresponds  to  a  thermalization  frequency  of  ~  2.4  x  10®s“^). 

Let  us  analyze  now  how  low-thermalization  of  beam  and  primary  electrons  affect  in  a  dif¬ 
ferent  way  the  plasma  response.  For  a  -C  1,  the  effective  secondary  yields  are 


a 

Ibp  —  _  r  72p) 

1  Osr 


l2p  ^ 


2Ti 

Eh  —  2T2 


(B.41) 


Therefore,  independently  of  the  emission  of  beam  electrons  (72^)  being  small  or  large,  low 
thermalization  makes  the  two  counter-streaming  SEE  beams  of  almost  equal  current  and  leads 
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Figure  B.6;  SEE  model.  Effect  of  Ssr{E)  =  6oexp{—E/Er)  for  6o  =  0,  0.2,  and  0.6,  Er/Ti  =  2, 
(1  —  SQ)Eh/Ti  =  2,  and  T2/T1  =  0.01.  Dashed  lines  correspond  to  use  the  approximate  function 
^sr{Ti)  of  Eq.(B.40).  Asterisks  correspond  to  the  CSL. 


to  an  effective  beam-to-primary  flux  ratio,  jbp,  very  small  [35].  The  consequences  are  Qp  ~  Qi 
and  qew  —  that  is  energy  losses  at  their  minimum.  On  the  other  hand,  the  sheath 

potential  fall  satishes 


^WQ  -  In 


(B.42) 


Here,  the  factor  a  is  measuring  the  depletion  of  the  tail  of  primary  electrons  impinging  the 
wall;  cr  -C  1  indicates  that  a  lower  (j)wQ  is  needed  to  adjust  the  electron  flux  to  the  ion  flux.  In 
addition,  -C  1  makes  that  SEE  beams  do  not  affect  the  sheath  potential  fall. 

The  behavior  of  4>wq  for  very  low  thernialization  is  not  immediate  since  Eq.  (B.49)  indicates 
that  gi  0  goes  to  zero  when  ^wq  0.  From  Eqs.  (B.33),  and  (B.47)-(B.49),  the  asymptotic 
plasma  behavior  for  the  ’very-low  thermalization  limit’(VLTL), 


a'  =  a{mi/meY^‘^  -C  1, 


(B.43) 


is 


neQ  ^  flpQ  ~ 


TrpQ  —  g  . 


(B.44) 


Figure  B.7  shows  the  evolution  of  these  parameters  in  the  low  and  very-low  a  ranges.  Notice 
that  as  long  as  a  >  ~  2  ■  10“^,  the  order  of  magnitude  of  the  main  plasma 

parameters  remains  unchanged.  Only  in  the  VLTL,  one  can  say  that  there  is  a  strong  depletion 
of  the  EVDF,  leading  to  T^pq/Ti  clearly  below  1.  It  is  interesting  to  observe  that,  for  the  VLTL, 
the  ion  velocity  at  the  entrance  of  the  sheath  satishes 


^i^riQ  —  ^'^rpQ- 


(B.45) 


Therefore,  at  the  VLTL,  the  Bohm  condition  is  preserved  but  the  relevant  temperature  is  the 
parallel  to  the  magnetic  held.  The  factor  3  is  due  to  the  strong  distortion  of  the  EVDF  shape 
from  a  Maxwellian  one  and  would  indicate  a  one-dimensional  adiabatic  behavior. 
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Figure  B. 7;  Low  thermalization  limit.  Same  four  cases  than  in  Fig.B. 3.  Dashed  lines  correspond 
to  the  very-low  thermalization  limit,  Eqs.  (B.44). 


B.6  The  charge  saturation  limit 


This  limit  is  the  transition  to  a  non-monotonic  structure  of  the  Debye  sheath  [54].  Its  impor¬ 
tance  here  lies  in  that  it  sets  a  local  minimum  of  the  sheath  potential  fall  and  a  maximum  of 
deposition  of  (dimensionless)  electron  energy.  The  CSL  is  dehned  by  condition  (B.50).  The 
numerical  computations  conhrm  that,  in  order  to  reach  the  CSL,  the  flux  of  primary  electrons 
must  be  much  larger  than  the  ion  flux  (i.e.  Qp/gi  3>  1),  which  means  that  the  beam  and 
primary  fluxes  are  very  similar,  that  is  jbp  —  1  in  Eq.(B.32).  The  ratio  Qp/ gi  is  plotted  in 
Fig.B. 8,  together  with  main  plasma  magnitudes  at  the  CSL;  an  asterisk  is  used  as  superscript 
for  parameters  at  the  CSL. 

The  temperature  of  primary  electrons  (Ti)  leading  to  sheath  charge-saturation  can  be  ob¬ 
tained  from  the  approximation  jbp  =  1-  Solving  this  equation  for  Ti  and  using  Eqs.  (B.  13)  and 
(B.40),  one  has 


2^1  ^  1  +  _  (i-^)(i-^:.)oT2 

El  a  I  -  o-  El' 


(B.46) 


with  6*^  =  Ssr{Ti).  This  expression  is  explicit  for  T^  only  if  the  yield  for  backscattered  electrons, 
6sr{Te),  is  constant.  Figure  B.8  and  Eq.(B.46)  show  that  the  plasma  temperature  T^  required 
to  reach  the  CSL  decreases  as:  a)  electron  thermalization  in  the  bulk  of  the  plasma  (i.e.  a) 
increases;  b)  the  average  SEE  yield  is  higher  (i.e.  Eb  is  lower);  or  c)  the  ’average’  temperature 
of  the  SEE  increases  (i.e.  either  6sr  or  T2  increase).  For  typical  values  of  Eb,  T2  and  Ti  in  Hall 
thrusters,  the  CSL  is  unattainable  for  low  thermalization  [say  for  6sr  >  0.3  and  a  <  0.20  (i.e 
\her/h  >  4.5)].  A  second  positive  aspect  of  low  thermalization  and  CSL  is  that  energy  losses 
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Figure  B.8;  Charge  saturation  limit.  Variation  of  main  parameters  with  a  for  6sr  =  0,  0.3,  0.6; 
other  parameters:  T2/T1  =  0.1.  Dashed  lines  for  correspond  to  Eq.(B.46). 


are  typically  one  order  of  magnitude  less  than  for  total  thermalization.  The  CSL  behavior 
obtained  here  is  the  same  obtained  in  Ref.  [35]  (for  6sr  =  0)  except  for  the  value  of  potential 
fall,  which  was  larger  there  since  the  depletion  of  primary  population  was  not  taken  into 

account. 


The  approximation  7^,^  =  1  cannot  be  used,  of  course,  to  compute  the  large  values  of 
parameters  (gp/gi)*  and  {qew /‘^Tigi)* ,  Eq.(B.32),  and  the  potential  fall  at  the  CSL, 
Eq.(B.33).  The  importance  of  computing  exactly  enough  {gp/gi)*  is  illustrated  in  Table  I, 
which  compares  the  results  of  the  Hobbs- Wesson  model  with  ours.  The  case  treated  by  these 
authors  corresponds  to  a  =  1,  6sr  =  0,  and  T2/T1  =  0,  except  that  they  computed  the  density 
of  primary  electrons,  Up,  using  the  complete  Maxwellian  distribution  (including  the  empty  tail 
of  electrons  collected  by  the  wall).  This  inconsistency  would  not  be  important  if  the  potential 
fall  at  the  CSL  was  large,  but  Table  I  shows  that  it  leads  to  errors  of  20%  on  {gp/gi)*  and 
the  deposition  of  electron  energy.  The  influence  of  the  temperature  T2  of  beam  electrons  on 
{gp/gi)*  is  also  worth  to  comment.  In  most  of  the  Hall  thruster  chamber,  the  ratio  T2/T1  is  of 
the  order  of  0.1,  which  makes  plausible  to  treat  the  ’beam  secondary’  population  as  cold  [9,  10] 
(except  in  the  very  vicinity  of  the  wall).  Errors  on  energy  deposition  and  sheath  potential  fall 
are  expected  to  be  of  order  T2/T1.  This  turns  out  to  be  true  except  for  the  CSL,  when  they  are 
of  order  ^/TfjT/^  as  Eq.(B.52)  explains  and  Table  I  illustrates.  Therefore,  ’small’  contributions 
(or  details)  of  the  electron  model  have  a  signihcant  effect  on  the  plasma  parameters  at  the  CSL. 
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{dp/ 9i) 

{deW /‘^TiQi)* 

^WQ 

9i 

FM,  T2 

=  0 

59.7 

59.7 

1.02 

1.08 

CM,  T2 

=  0 

74.1 

74.1 

0.91 

1.03 

CM,  fa 

=  0.01 

79.8 

79.0 

0.82 

1.04 

CM,  fa 

=  0.1 

94.5 

85.1 

0.63 

1.06 

TABLE  1.  CSL  parameters.  FM  stands  for  the  full-Maxwellian  distribution  assumed  by 
Hobbs  and  Wesson;  CM  stands  for  the  corrected-Maxwellian  with  the  empty  tail  of  collected 
electrons.  Other  parameters:  cr  =  1  and  6sr  =  0. 


Ti  (eV) 


Figure  B.9:  Main  plasma-wall  parameters  (sheath  potential  fall,  ion  current  to  wall,  and  electron 
energy  deposition  at  walls  )  as  function  of  the  perpendicular  temperature  of  primary  electrons 
Ti  for;  different  thermalization  levels  (cxi  =0.01,0.1,0.3,  and  1);  n^Q  =  3  ■  10^^m“^;  Ei  =  40eV; 
T2  =  2eV;  and  6sr  =  0  (solid  lines),  =  0.45  and  Er  =  50eV  (dashed  lines).  The  region 
to  the  right  of  the  asterisk  for  the  solid  curve  of  a  =  1  corresponds  (approximately)  to  the 
charge-saturation  regime. 
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B.7  A  practical  example  and  additional  considerations 

Figures  B.3  to  B.8  have  presented  dimensionless  curves  which  allow  to  evaluate  the  relevance  of 
the  different  aspects  of  the  model.  In  order  to  complete  the  illustration  of  the  influence  of  partial 
thermalization  on  the  plasma-wall  interaction,  Fig.B.9  presents  dimensional  results  for  typical 
plasma  values  of  a  conventional  stationary  plasma  thruster(SPT)  of  the  l-2kW  range  [63]  and 
lateral  walls  made  of  Boron  Nitride;  the  SEE  data  for  this  material  has  been  taken  from  Ref.  [55]. 
A  crossover  energy  Ei  of  40  eV  is  taken,  and  solid  and  dashed  lines  differ  on  the  contributions 
of  true-secondary  and  elastically-reflected  electrons.  Thermalization  levels  from  a  =  0.01  to 
a  =  1  are  considered;  for  a  chamber  width  of  h  =  15mm  and  electron  energies  of  50eV,  the  range 
a  =  10“^  —  10“^  corresponds  to  mean  free  paths  of  0.14  —  1.4m  and  thermalization  frequencies 
of  2  ■  10®  —  2  ■  10^s“^.  Figure  B.9  plots  the  sheath  potential  fall,  the  ion/electron  current  to 
the  wall,  and  the  electron  energy  deposited  at  the  wall;  ion  energy  losses  follow  the  trends  of 
the  sheath  potential  fall.  Notice  that  particle  and  energy  fluxes  to  the  wall  are  proportional  to 
the  plasma  density  at  the  sheath  edge,  UeQ  (which  is  about  a  50-60%  of  the  electron  density 
at  the  channel  median).  The  plots  shows  that  the  thermalization  level  is  very  crucial  when 
estimating  plasma-wall  parameters.  In  particular  electron  energy  losses  can  vary  by  two  orders 
of  magnitude,  but  for  the  range  of  temperatures  of  interest,  the  large  electron  energy  losses  at 
the  charge-saturation  regime(CSR)  are  attained  only  if  the  thermalization  level  is  high.  On  the 
other  hand,  the  presence  of  a  signihcant  fraction  of  ’hot’  backscattered  electrons  (i.e.  dashed 
lines)  has  a  secondary  effect,  except  for  the  retardation  of  the  CSL  at  high  thermalization. 
These  are  positive  news  for  modelling  since  the  uncertainties  on  the  distribution  function  of 
secondary  electrons  are  large.  The  values  of  ion  current  to  the  wall  in  Fig.B.9  are  of  the  order 
of  those  measured  by  Kim  et  ah  [63].  As  explained  before,  the  ion  current  into  sheath  and 
wall  depends  mainly  on  the  plasma  physics  in  the  bulk  region,  and  is  weakly  affected  by  the 
thermalization  level  and  the  SEE  details.  The  reduction  of  the  ion  current  for  a  =  0.01  is  due 
to  the  a  decrease  of  the  effective  parallel  temperature  of  primary  electrons;  the  increase  of  the 
ion  current  for  a  =  1  at  the  CSR  is  the  effect  of  an  increase  of  the  secondary  beam  density  (see 
UbQ  in  Fig.B.3). 

The  curves  for  <7  =  1  and  6sr  =  0  in  Fig.B.9  correspond  to  the  model  presented  in  Ref. [10] 
and  implemented  now  in  the  simulation  code  HPHall-2  of  the  plasma  discharge[3].  There  are 
no  difficulties  on  implementing  the  hve-parameter  SEE  model  proposed  here  in  the  simulation 
code,  but  we  can  predict  that  changes  are  going  to  be  limited.  However,  the  implementation 
of  a  partial  thermalization  model  for  electrons  in  the  simulation  code  is  far  more  challenging. 
The  electron  fluid  model  currently  implemented  in  HPHall-2  (and  in  other  hybrid  codes)  for  the 
bulk,  quasineutral  region  is  based  on  a  one-temperature,  quasi-Maxwellian  electron  distribution. 
An  electron  fluid-type  model  for  low  thermalization  requires  to  assume  (i)  a  two-temperature 
distribution  for  primary  electrons,  and  (ii)  an  independent  distribution  for  the  true-secondary 
electron  beams.  This  implies,  at  least,  to  postulate  (i)  conditions  relating  the  parallel  and 
perpendicular  temperatures  of  primary  electrons,  (ii)  equations  for  the  macroscopic  magnitudes 
of  beam  electrons,  and  (hi)  source  terms  ’transforming’  beam  electrons  into  primary  ones.  An 
additional,  important  challenge  is  that  the  fluid  equations  for  magnetized  electrons  must  be 
expressed  in  the  non-cartesian  reference  frame  linked  to  the  curved  magnetic  streamlines. 

Therefore,  we  believe  that  there  are  several  important  issues  related  to  the  evolution  of 
secondary  electrons  that  remain  largely  unsolved.  These  should  be  understood  prior  to  affording 
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the  implementation  of  a  new  electron  model  in  codes  such  as  HPHall-2.  First,  the  determination 
and  evaluation  of  the  main  physical  processes  that  thermalize  (or  isotropize)  the  EVDF  are 
poorly  known.  Second,  Sydorenko  et  ah  [64]  have  shown  recently  that  the  two-stream  instability 
can  ’heat’  and  trap  an  important  fraction  of  the  secondary  beams,  thus  indicating  the  existence 
of  a  trapping  process  for  secondary  electrons  independent  of  the  thermalization  of  the  primary 
EVDF  (modelled  through  a).  Third,  a  consistent  model  for  trapped  secondary  electrons  and 
their  ’thermalization’  into  primary  electrons  (in  a  slower  time-scale)  must  be  derived.  And 
fourth,  the  fluid  model  of  the  secondary  beams  should  take  into  consideration  the  different 
magnetic  effects  acting  on  them,  such  as  the  E  x  B  and  curvature  drifts,  magnetic  mirror 
trapping,  and  the  oblique  incidence  of  the  magnetic  held. 


B.8  Conclusions 

A  model  of  the  plasma  interaction  with  the  dielectric  walls  of  a  Hall  thruster  has  been  presented. 
It  accounts  for  partial  thermalization  of  the  electron  population  through  a  single  parameter 
(7  and  includes  a  two-population,  four-parameter  model  for  SEE.  Analytical  expressions  are 
obtained  for  the  main  parameters  characterizing  that  interaction,  such  as  the  particle  and 
energy  huxes  to  the  walls  and  sheaths,  which  are  needed  as  boundary  conditions  of  quasineutral 
models  of  the  full  discharge. 

The  behavior  for  low  thermalization  is  shown  to  differ  greatly  from  the  commonly-used, 
high-thermalization  case.  This  is  very  relevant  for  Hall  thrusters,  where  there  is  a  growing 
conviction  that  electron  thermalization  is  low  (at  least  for  primary  electrons).  At  low  thermal¬ 
ization,  energy  losses  are  close  to  its  minimum,  the  charge  saturation  limit  is  not  attainable, 
and  the  sheath  potential  is  small;  the  different  roles  of  beam  and  primary  electrons  on  these 
characteristics  have  been  analyzed.  Signihcant  decreases  of  the  parallel  temperature  of  primary 
electrons  and,  therefore,  of  the  plasma  flux  into  the  sheath  (through  fulhlment  of  the  Bohm 
condition)  take  place  only  at  the  very-low  thermalization  limit  {a 

The  investigation  of  the  emission  model  for  secondary  electrons  has  shown  that  the  presence 
of  a  relevant  fraction  of  elastically  reflected  electrons  affects  greatly  the  response.  They  tend 
to  amplify  the  relative  densities  of  untrapped  electrons;  their  effect  on  the  net  primary  and 
beam  fluxes  comes  out  from  the  zero  electrical  current  balance.  It  is  reiteratively  found  that 
the  role  of  the  sheath  potential  fall  is  to  adjust  the  primary  electron  flux  to  wall  and  not 
viceversa.  Although  most  of  the  analysis  is  carried  out  for  an  energy-independent  yield  of 
reflected  electrons,  a  temperature-dependent  yield  expression  is  proposed,  which  avoids  integrals 
expressions  at  the  same  time  that  it  recovers  approximately  the  reduction  of  that  yield  with 
the  impact  energy. 

The  charge  saturation  limit  is  attained  when  net  beam-to-primary  net  fluxes  ratio,  'jbp,  is 
very  close  to  one,  that  is  the  ratio  between  the  primary-electron  to  ion  flux  ratio,  Qp/pi  = 
(1  —  7f,p)“^,  is  very  large.  The  CSL  requires  larger  plasma  temperatures  as  the  thermalization 
decreases  and  the  average  yield  of  reflected  electrons  increases.  Additionally,  the  sensitivity  of 
the  value  of  Pp/ Pi  at  the  CSL  to  small  variations  on  the  model  properties  has  been  stood  out. 
As  examples  we  have  shown  that  the  influence  of  a  small  emission  energy  of  true-secondary 
electrons  (T2)  or  the  (inconsistent)  inclusion  of  the  collected  tail  of  primary  electrons  are  less 
marginal  than  what  could  be  expected. 


no 
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Appendix:  electron  magnitudes  in  the  sheath 


The  electron  density  ne{4>)  =  f  dvf{(j),v)  has  contributions  from  primary  and  beam  electrons. 
Using  Eq.(B.20)  for  the  EVDF  and  writing  he  =  hp  +  hf,,  one  has 


n. 


,(0)  =  e'^ 


(1  +  Ssr)o'n 


erfc\/  0  —  (pw 


Tb  fn. 


nb{4>)  =  YyY^  -  (t)w)/T2, 


(B.47) 


with  /ift  =  (2  -  (t)(1  -  5sr)nibp- 

The  parallel-to-perpendicular  temperature  ratio  of  the  primary  population  at  point  Q,  is  of 
interest  also.  It  is  readily  seen  that  the  perpendicular  temperature  is  Ti.  Then,  the  temperature 
ratio,  TrpQi  satishes 
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mvlfpQ 
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d^v 


1  - 


(i  +  Uy-i 
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Pj-1/2^  +  eriC(/>^Q 
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where  the  integration  domain  is  the  whole  velocity  space  and  the  small  contribution  of  the 
macroscopic  radial  kinetic  energy  of  electrons  has  been  neglected. 

Using  the  above  expressions  for  the  electron  densities,  the  sonic  Bohm  condition  (B.37)  for 
gi  becomes 
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Finally,  the  3  contributions  to  the  CSL  condition. 


hfe]  dc() 


0, 


(B.50) 


are 


j^fiidlj)  =  [Jg^f  +  ‘i^wQnlQ-  go 

V  /  “eQ 

jw  =  ^pQ  ~  ^pw  -  - — (B.51) 
J^fibd^  =  {ubQ  -  nbw)^^  + 


For  T2/T1  -C  1  (and  (pwQ  ~  1);  one  has 
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